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Preface 


The need for increased efficiency in the use of our 
energy resources has stimulated applied research in many 
areas. Recently progress has been made in the field of 
aerodynamics, where the development of the supercritical 
wing promises significant savings in consumption 

of aircraft operating near the speed of sound. Computa~ 
tional transonic aerodynamics has proved to be a useful 
tool in the design and evaluation of these wings. 

We present here a numerical technique for the design 
of two-dimensional supercritical wing sections with low 
wave drag. The method is actually a design mode of the 
analysis code r developed by Bauer, Garabedian, and Korn 
[2,3,4]. This analysis code gives excellent agreement 
with experimental results and is used widely by the air- 
craft industry. We hope the addition of a conceptually 
simple design version will make this code even more useful 
to the engineering public. 
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I . INTRODUCTION 


Description of the Problem 

In this section we discuss seme of the principles 
behind the supercritical wing and describe our contribu- 
tion to the subject. 

General considerations show that the range of an 
aircraft is roughly proportional to the parameter M L/D. 
where L is the lift, D is the drag, and the free stream 
Mach number is the ratio of the aircraft's speed to 
the speed of sound. The top curve in Figure 1 shows the 
general behavior of thxs parameter as the Mach number M 

CO 

is varied. The value of M^L/D that maximizes the range 
aircraft occurs near a region of rapid increase 
in drag known as drag rise, shown by the bottom curve of 
Figure 1. At this speed the flow is observed to be 
transonic, with regions of supersonic flow appearing where 
the air accelerates over the wing, when the free stream 
Mach number has the value depicted in Figure 1, the 
maximum speed of the flow is equal to the speed of 
sound. When > .M^ the flow is said to be supercritical. 

Figure 2 illustrates some important characteristics 
of the flow past an airfoil at speeds corresponding to drag 
rise. The region of supersonic flow is terminated by a 
shock, where the pressure is observed to bo discontinuous. 
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As the speed of the wing is increased, the supersonic zone 
grows in size and the pressure jump becomes larger. The 
occurrence of such shocks in the flow imposes a retarding 
force on the wing known as wave drag, which is one reason 
for the drag rise seen in Figure 1. The large pressure 
gradients present in strong shocks can also induce separa- 
tion of the boundary layer of air that adheres to the wing 
because of friction; separation results in a decrease in 
lift and more drag. The study of transonic flow is 
therefore important not only because transonic flow 
encompasses the most economical regime for aircraft opera- 
tion, but because the deterioration of an aircraft's 
efficiency at higher speeds is due to transonic effects. 

The supercritical wing is designed to delay the onset 
of drag rise to higher Mach numbers . Since the efficiency 
of the wing is governed by the optimal value of M^L/D, 
postponing the onset of drag rise to higher Mach numbers 
results in a corresponding decrease in the fuel requirements 
of the aircraft. The delay in drag rise can be effected by 
constructing the wing so that the strong shocks accompanying 
supersonic zones of moderate size on conventional wings are 
replaced by weaker shocks with less wave drag and no appreci- 
able boundary layer separation. 

We are mainly concerned here with the contributions to 
supercritical wing technology made by computational transonic 
aerodynamics. The numerical solution of the partial 
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differential equations of gas dynamics provides a theoreti- 
cal means for both the design and evaluation of supercritical 
wings. For example, two-dimensional shockless airfoils can 
be obtained by calculating real analytic solutions to the 
hodograph equations of transonic flow. These airfoils have 
the property that at a specified speed and angle of attack, 
the calculated two-dimensional transonic flow is smooth. 

This guarantees that the wave drag will be small near at 
least one operating condition. It may happen that there is 
an increase in wave drag at supercritical Mach numbers below 
the design condition known as drag creep. Since a wing must 
operate efficiently over a range of conditions it is desir- 
able to avoid the occurrence of drag creep. Provided this 
is done, a practical approach to the supercritical wing is 
to design the wing using computer-generated shockless airfoils 
in each cross-section. 

It is also possible to evaluate the performance of Vv’ings 
at off-design conditions using computer codes. Programs 
that calculate the three-dimensional transonic flow past 
wing-body combinations are used regularly by the aircraft 
industry. Considerable savings can be realized by replacing 
preliminary wind tunnel testing of new wing designs by such 
computer simulation . 

There is presently much interest in the possibility of 
developing a numerical scheme fot the design of three- 
dimensional wing-body combinations. Hodograph methods employed 
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in the design of two-dimensional shockless wing sections 
cannot be used for this purpose. Our contribution in 
this direction is a two-dimensional design code based on 
techniques that may prove useful in such an endeavor. 
Although our method does not produce shockless airfoils, 
we show that it is possible to obtain wing sections with 
low wave drag by using an artificial viscosity to smear 
shocks properly. The design procedure is outlined in 
Section 4.1. 

Our program is actually a new design mode of the 
two-dimensional analysis code H developed by Bauer, 
Garabedian, and Korn [2,3,4]. This analysis code solves 
the direct problem of obtaining the transonic flow past 
a given wing section. The code includes a turbulent 
boundary layer correction which gives a reliable approxima- 
tion to the drag due to skin friction and predicts boundary 
layer separation. Drag estimates obtained with the 
code are in good agreement with experiment, and the program 
has found wide acceptance in the aircraft industry. 

There are two major steps in the operation of the 
analysis routine. Given the coordinates of the airfoil, 
the region exterior to the airfoil is mapped conformally 
to the interior of the unit circle as in Sell's treatment 
of suocritical flow past 'in airfoil [33] . The nonlinear 
partial differential equ *tions of transonic flow are then 
solved iteratively in the unit circle using a type- 
dependent difference scheme similar to the one first 
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de /eloped by Murman and Cole [27]. If boundary layer 
corrections are desired, the shock wave - boundary layer 
interaction is simulated by iterating between inviscid 
flow calculations and boundary layer corrections until 
convergence is achieved. 

The design modification we have added to the code 
solves the inverse problem of calculating the shape 
an airfoil must have in order to achieve a specified 
pressure distribution. When operating in this mode, an 
initial guess is provided for the shape of the desired 
airfoil and the region exterior to this airfoil is mapped 
into the unit circle as in the analysis mode. A number 
of flow iterations are performed using an artificial 
viscosity that inhibits the formation of shocks, as 
described in Section 2.2. The pressure disti.ibution 
resulting from these calculations is then compared with 
the desired input pressure distribution and a better 
approximation to the desired airfoil is obtained, as 
described in Section 2.3. This new profile is mapped 
to the unit circle as in the first step and the process 
is repeated until the approximations converge. A boundary 
layer correction may then be calculated on the basis of 
the last pressure distribution. The desired airfoil is 
obtained by subtracting the displacement thickness of the 
turbulent boundary layer from the coordinates of the 
computed profile. 
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The inverse method transfers the difficulty in 
designing wings from determining the coordinates of the 
airfoil to finding pressure distributions that give rise 
to airfoils with desired specifications. We therefore 
include descriptions of some pressure distributions that 
generate airfoils with low wave drag, and indicate how 
to modify the input distribution in order to obtain air- 
foils with a desired lift, thickness-to-chord ratio, and 
design Mach number. 

The remainder of the paper is organized as follows. 

The mathematical statement of the problem is formulated 
in Chapter II . The computational procedure is outlined 
in Chapter III. In Chapter IV we present results obtained 
with the design mode, together with some comparisons to 
airfoils obtained by other methods . Chapter V is more 
theoretical in nature and includes a convergence proof for 
the design problem in a special case of subsonic flow. 

We provide a description of the modified version 
of the Bauer, Garabedian, and Korn analysis code H in 
Chapter VI. 

I would like to express my gratitude for the advice 
and encouragement of Paul R. Garabedian, who suggested this 
problem and made the work possible. I am also grateful for 
the help of Frances Bauer and Antony Jameson at various 
stages of the research, and for a fast and accurate typing 
job by Connie Engle. 
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2* References to Other Wnri» 


Transonic flow research has a colorful history [5,29]. 
In the late 1940's. arguments to the effect that smooth 
transonic flows past arbitrary profiles should not generally 
be expected to exist were formulated by Buscmann [lOJ , 

Frankl [I 5 ], and Guderley [16]. These observations raised 
doubts about the physical significance of the smooth solu- 
tions to the steady, two-dimensional potential equation for 
transonic flow that were known at the time, it xvas observed 
experimentally that transonic flows generally exhibit shocks 
when the supersonic zones are of moderate size, but there 
were occasional instances of near-shockless flow that 
seemed to contradict the implications of the nonexistence 
theorems. A "transonic controversy" developed over the 

true nature of transonic flows in general and of shockless 
flews in particular. 

The controversy attracted considerable attention from 
mathematicians in the hopes that a rigorous investigation of 
whether the flow problem was well-posed would help clarify 
matters. From this viewpoint, a satisfactory demonstration 
that the problem of finding smooth transonic flows past 
convex symmetric profiles was not correctly set was supplied 
by Morawe^z [ 26 ], although the apparent discrepancies between 
theory and experiment remained unresolved. 

More progress was made in the early I960's with the 
experimental work of Pearcey f 31] , v;ho was able to systemati- 



cally produce near-shockless flows past wing sections having 
a suction pressure peak near the nose of the profile. The 
subsequent development of numerical techniques capable of 
treating transonic flows with shocks brought about a 
reinterpretation of the nonexistence theorem‘s since the 
computational problem seems to be correctly 
set in terms of weak solutions. Results of both 
experiment and computation show that snockless and 
neighboring near-shockless solutions do in fact have 
physical significance and can provide an important means 
of reducing the drag experienced by airfoils travelling 
at transonic speeds. 

Several numerical techniques have been developed for 
the design of two-dimensional supercritical wing sections 
using inviscid flow theory. We can distinguish between 
approaches relying on hodograph methods and the remaining 
approaches . 

The hodograph transformation conrists of reversing 
the roles of the dependent and independent variables in the 
flow equations with the result that the partial differential 
equations are linear. Using this transformation, several 
methods have been devised to allow the systematic computa- 
tion of airfoils that have shockless flows at given operat- 
ing conditions [2,4,8,30]. Such airfoils necessarily have 
low wave drag at nearby operating conditions, although 
drag creep can occur elsewhere. 



other approaches to the design problem do not 
generally provide shock-free solutions to the equations. 
This is not necessarily a disadvantage, since for some 
applications it is possible that an airfoil designed with 
a weak shock might have an overall performance that is 
as good or better than a shockless airfoil with similar 
specifications. The main difficulty is to find pressure 
distributions that will generate airfoils with low drag 
levels . 

Some design methods use the approximations of small 
disturbance theory and thin airfoil theory [11,22,31]. 

With this approach the solution is expanded in terms of a 
parameter describing the thickness of the profile, which 
is assumed to be small. This has the advantage that to 
leading order the profile can be replaced by a given slit. 
The desired pressure distrubution along the surface of 
the airfoil can then be used in a boundary condition applied 
at the slit, and so the difficulties caused by the unknown 
boundary are avoided. The coordinates of the desired air- 
foil can be determined from the resulting velocity components 
This technique has the disadvantage that the flow is not 
represented correctly near the stagnation point at the 
leading edge of a Munt-nosed airfoil. it should be noted 
that applications of this technique to the design of three- 
dimensional wings have been initiated [17,34]. 





The design methods of Carlson tl2] and Tranen [36] 
solve the inverse problem for the full potential equation 
with a free boundary. Both of these methods use the 
prescribed pressure dxstribution in a boundary condition 
for the determination of the velocity potential, and 
calculate the position of the surface of the profile by 
using the condition of flow tangency along the body. In 
Carlson's method the calculation is performed using Cartesian 
coordinates. The coordinates near the nose are given in 
advance and the remainder of the profile is determined as a 
free boundary. Tranen uses the analysis code H to perform 
the flow calculations and to provide a computational domain, 
and proceeds by alternating between analysis and design 
computations. At each c^cle the user modifies the prescribed 
pressure distribution in order to achieve convergence. 

Ihe design procedure of Hicks and his associates [18] 
is based on the use of a numerical optimization routine 
together with the analysis code H to minimize the drag coeffi- 
cient with respect to design variables that describe 
the shape of the profile, while satisfying various constraints 
on the operating conditions and geometry. This technique 
has the advantage of drag reduction without the necessity 
of choosing the pressure distribution. Its main drawback 
IS the large amount of computing time required when many 
parameters are allowed to vary. 
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The method we present for supercritical wing design 
uses the prescribed pressure distribution in a boundary 
condition for the determination of the conformal mapping 
from the unit circle to the desired airfoil. This 
approach to the desi;a problem is similar in spirit to 
Lighthill's inverse method [23], which is based on the 
linear theory of incompressible flow and so does not 
require an iterative procedure to determine the flow and 
profile. Other incompressible treatments along these 
lines have also appeared [ 1 , 13 ] . 

The practical success of an inverse method of airfoil 
design depends on the prescribed pressure distribution. 

It is therefore important to study the relation between 
the assigned pressure distribution and the performance of 
the resulting airfoil [7,29]. Much work remains to be done 
on this aspect of the problem. The many shockless flows 
produced by hodograph methods provide a good basis for 
investigation . 
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II. THE PARTIAL DIFFERENTIAL EQUATIONS OF TRANSONIC FLOW 


In this chapter we consider the mathematical formula- 
tion of the design problem. We summarize the basic equa- 
tions of motion for gas dynamics and discuss the boundary 
conditions appropriate for the direct and inverse problems. 

1 . The Equations of Gas Dynamics 

We begin /ith some comments about the choice of equa- 
tions to describe the problem. We wish to model the 
flight of aerodynamically efficient wings at transonic 
speeds. We are especially concerned with the drag on such 
bodies, which includes forces due to skin friction and 
shocks. For our treatment to be of practical use we must 
consider equations which allow estimates of the wave drag 
due to shocks, and we must provide for the calculation 
of viscous effects. Furthermore, we must choose equations 
that are compatible with the inherent storage limitations 
of computers . 

It is common in experiment as well as theory 

to treat the case of steady, two-dimensional flow past a 
wing of uniform cross section. This provides a good 
approximation of the flow near the middle section of 
a three-dimensional wing with a straight leading edge 
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moving with a constant velocity. This geometrical simpli- 
fication permits the use of just two independent var;=»bles, 
which we take to be the x and y coordinates. 

Another important simplification is possible if the 
sii'foil is streamlined so that viscous effects are confined 
to the immediate vicinity of the profile. In this case 
the flow outside the boundary layer can be obtained from 
lower order partial differential equations describing 
inviscid fluid motion. The inviscad solution can then 
be used to calculate a boundary layer correction to the 
flow past the airfoil, making it possible to obtain esti- 
mates of the drag due to skin friction [28,32]. Separation 
of the boundary layer should be avoided for aerodynamical 
reasons, too, so it is important for the theory to give 
reliable estimates of the growth of the boundary layer. 

Inviscid fluid flow can be described by conservation 

laws consisting of nonlinear, first order ’partial differ- 

; 

ential equations involving the velocity components of the 
flow and two thermodynamic variables such as the density 
and entropy. The conservation laws also provide shock 
conditions which determine the jump in these quantities 
across a surface of discontinuity in the flow. For the 
case of flows past thin bodies at speeds close to the 
speed of sound, the shocks are usually weak in the sense 
that che ^ump in velocity across the shock is small compared 
to the soeed of sound. The jump in entropy across a shock 
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is of third order in the shock strength; to a good 
approximation, the change in entropy can therefore be 
neglected in a weak shock. With this assvimption, the 
entropy of a fluid particle is constant throughout its 
motion, and the flow may be considered isentropic. 

As a result of considering the entropy to be con- 
served across a shock, the shock condition expressing 
conservation of the normal component of momentum is lost. 

The defect in this quantity can be interpreted as an 
approximation of the wave drag exerted on the airfoil 
by the shock. 

For isentropic flow, the velocity field is irrotational 
if the flow is uniform at infinity. This permits the 
introduction of a velocity potential whose derivatives are 
the velocity components. The inviscid equations of motion 
can then be reduced to a single second order partial differ- 
ential equation for the potential, and in the computation 
it is only necessary to store values of a single dependent 
variable. 

We shall list below the equations describing the flow 
of an inviscid, isentropic gas. In Section 3.4 the equa- 
tions used to describe a turbulent boundary layer correction 
will be discussed. It is found in practice that these 
equations provide a description of the transonic flow past 
an airfoil that agrees well with experiment over a wide 
range of conditions [3,4]. 
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The equations describing the steady two-dimensional 
motion of an ideal polytropic gas are familiar [14,25], 

From thermodynamics we have the equation of state 

(2.1) p = A(s)p^ 

where p, p, and S are the pressure, density, and specific 
entropy of the gas. A(s) is a known function of the 
entropy and y > 1 is a constant depending on the nature 
of the gas. 

Conservation of mass gives 

(2.2) (pu) + (pv) = 0 

A y 

where u and v are the x and y components of the velocity. 
Similarly, the conservation of momentum asserts that 


(2.3) 

(UU^ + VUy) 

+ P. 

(2.4) 

(UVj^ + Wy) 

+ P; 


and the conservation of energy gives 

(2.5) uS^ + VS^^ = 0 . 

X y 

As mentioned above, we consider the case of constant 
entropy, so that (2.5) is automatically satisfied and A(S) 
in (2.1) is a constant. The flows we consider become 
uniform at large distances. It then follows from a theorem 
of Kelvin that the flov/ is irrotational. 
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( 2 . 6 ) 


= 0 . 


U - V 

y X 


Using (2.1), (2.2) and (2.6) we obtain Bernoulli's law 


(2.7) 


«2 + v2 


- c* , 


Y-1 2 Y - 1 


where c - dp/dp is the square of the local speed of 
sound and c* is a constant known as the critical speed. 
The dimensionless ratio 




is called the local Mach number and is greater than one if 
2 2 _ 2 2 

u + V q > c* (locally supersonic flow) and less than 
one if q < c* (locally subsonic flow) . 

According to (2.2) and (2.6), there are two functions 
and iji such that 


(2.8a) 

(2.8b) 


u 


= ’^y/P 


V = = - VP 


<P and are the velocity potential and stream function of 
the flow. We may obtain a single equation for « from (2.1), 
(2.2), and (2.7), 

(2.9) (c^ - + (c2 - = 0 . 

This is the partial differential equation that is the basis 
of our numerical work, ip satisfies a similar equation. 
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Equation (2.9) is a quasilinear partial differential 
equation which is elliptic when < i a„d hyperbolic when 
M >1. We are interested in the case of transonic flow, 
so that (2.9) has mixed type in the region of interest. 

It is appropriate to consider weak solutions to (2.9) 
or (2.2) under the conditions that ^ is continuous and 
that any shocks present are compressive. This corresponds 
to tne proper entropy inequality in nonisentropic flows. 
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2 


The Direct Problem 


In this section we discuss the formulation of the 
direct, or analysis, problem of determining the flow 
past a given wing section. 

We consider an airfoil with coordinates (x(s),y(s)j 
parametrized by arc length s measured from tail to tail 
as in Figure 2. The included angle at the trailing edge 
is denoted by e . The frame of reference is chosen 
so that the airfoil is at rest and the air has a resulting 
velocity u = (q^ cos a, sin a) at infinity, where ^ 
the angle of attack a gives the direction of motion relative 
to fixed coordinate axes. 

The fact that the flow must be tangential to the 
surface of the airfoil provides one boundary condition 
for (2.9). If V is a unit normal to the profile, then 
this condition can be stated as 

(2.10) 0-u = H = 0 

- u V 

on the curve (x(s),y(s)). Since the airfoil is a stream- 
line of the flow, this fact can also be expressed by 

(2.11) >1' (x (s) ,y (s) ) = constant. 

We consider lifting profiles with cusped trailing edges, 
0 £ E << 1, in which case the potential <|) need not be 
single-valued. The circulation T = [(J>] of the flow around 
the airfoil is then uniquely determined by the Kutta- 
Joukowski condition 
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( 2 . 12 ) 


|u(x(0) ,y(0) ) I < oo 


which states that the velocity must be finite at the 
trailing edge. „ s > o the flow necessarily nas a 
stagnation point at the trailing edge. This is not the 

0. The presence of a stagnation point at 
the tail of the airfoil should generally be avoided 
Since the resulting adverse pressure gradient promotes 
separation of the boundary layer. 

It can be shown ( 24 1 that t has an asymptotic 
expansion 

(2.13) t V cos( 0 -„) . tan-l(B tan(9-c)) 

as r^ = . y2 . ,, " 1 - a.nd 9 = tan'l y/r. 

This is similar to the corresponding er.iansion for 

incompressible flow. The Prandtl-Glauert scale factor B 

commonly occurs in linearized treatments of compressible 
flow. 

analysis problem it is convenient to express 
Bernoulli's law (2.7) in the form 

(Y-Dm;^ 

where is the Mach number at infinity. „tth this nota- 
tion, the direct problem can be formulated by specifying 

the airfoil coordinates (x(s),y(s)), the angle of attack n. 


(2.14) 


— y + si - _2fi I ^ 

2 ^ ^ 2 If 

'■ (Y-l)Mf J 
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and < 1. in this case we are free to choose the units 
so that is normalized to one. The flow is then obtained 
by solving the partial differential equation (2.9), along 
with the boundary conditions ( 2 . 10 ), ( 2 . 12 ), ( 2 . 13 ), and 

(2.14) . 

The lift of the airfoil is proportional to the 
circulation r. When the angle of attac)« as varied, the 
circulation of the flow adjusts so that the Kutta condi- 
tion (2.12) is satisfied. For a given Mach number M 

oo ' 

the lift is therefore a function of the angle of attacJc. 

A variant of the above formulation of the problem that 
is useful in applications is to specify the lift of the 
airfoil instead of a. The angle of attack necessary to 

pro .u .e this 3ift is then determined by imposing the Kutta 
condition. 
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3. 


The Inverse Problem 


In this section we describe the inverse, or desigi^ 
problem of calculating the shape that a profile must have 
in order to achieve a given pressure distribution. 

Suppose there is a flow past a profile such as the one 
depicted in Figure 2. If the velocity of the flow along 
the profile is written in the form 

u(s) - iv(s) = Q(s) , 

then the direct problem consists of specifying the angle 
0(s), which determines the shape of the airfoil, and solv- 
ing for the flow. For the inverse problem, the values of 
Q(s) are given and both the flow and the body are to be 
^sts^rnmed. Since Bernoulli's law (2.7) provides a corres- 
pondence between the values of and c^ = constant 
we may formulate the inverse prob.lem in terms of either 
*3 P' and the choice of q is only a matter of mathematical 
convenience. 

The inverse problem is seen to be a free boundary problem 
and IS for this reason more complicated than the direct 
probl€.m. The coordinates (x(s),y(s)) of the airfoil are 
now unknovm and are to be determined from the knowledge 
of the velocity distribution Q(s) . We may write the 
relationship between (}> and Q(s) as an additional boundary 
condition 
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t . (2.15) <f>(x(s) ,y (s)} = Q(s) 

f interval 0 < s < £, where I is to be the total 

I length of the airfoil. 

* It is also necessary to specify the constant in 

Bernoulli's law. in the direct problem the Mach number 

t 

I and speed at infinity are prescribed as in (2.14); for 

the design problem we give instead the value of the 
critical speed c* in (2.7). This means that the Mach 
numbers of the flow along the airfoil are specified. 

I The choice of c* determines the type of the equation (2.8) . 

If > max |Q(s)|, the flow will be subsonic and (2.9) 
will be elliptic; if there are points with |Q(s)| > c* , 
the flow will be transonic and (2.9) will have mixed type. 

We remark that with this formulation is not speci- 
fied as data, but must be determined along with and 
(x(s) ,y(s)) . 

The fact that Q(s) is to be the velocity distribution 
of a flow past an airfoil places restrictions on the form 
Q(s) may have. A typical choice for Q(s) is illustrated 
in Figure 3. Q(s) must have a zero corresponding to the 
stagnation point that forms at the nose of the airfoil, 
and then must be nonzero along the surface of the airfoil 
until the tail is reached. At tlie tail the velocity must 
be continuous and may be taken to be nonzero provided the 
included angle e at the trailing edge is zero. 


22 



Q(s) must satisfy further compatibility conditions in 
order to determine profiles defined by simple closed curves. 

Note that since the circulation of the flow is given 
by the integral of the speed along the profile, the lift 
of the airfoil can be calculated from the prescribed 
velocity distribution. The angle of attack may still be 
specified in the asymptotic form (2.13), since in the design 
problem the airfoil is free to rotate with respect to the 
coordinate system so as to satisfy (2,12) 

To summarize, the design problem is posed by specify- 
ing the speed distribution Q(s), the critical speed c* , 
and the angle of attack a. The potential of the flow and 
the airfoil coordinates are then obtained by solving the 
equations 


(2.9) 
(2.7) 

( 2 . 10 ) 

(2.15) 

( 2 . 12 ) 

(2.13) 


2<J>xVxy + “^y^^^yy 

d)^ + (J)^ 2 

!l + ^ 1 Y+l ^2 

2 Y - 1 2 ' 


= 0 , 


Iv (x(s),y(s)} = 0 ; 


^ <^(x(s) ,y(s)j = Q(s) ; 


|u(x(0) ,y(0)j I < oc 


4> q„r cos(e-a) + 


tan"^[6 tan(G-a)) 
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III. DISCUSSION OP THE COMPUTATIONAL PROCEDURE 

In this chapter wc describe the method used to solve 
the design problem outlined in Section 2.3. m Section 3.4 
we also provide a brief suitunary of the equations used to 
compute the boundary layer correction. 


Overvie w of the Computation 

The procedure we describe here is based on the Bauer, 
Garabedian, and Korn analysis code H [2,3,4] which solves 
the direct problem described in Section 2.2. The analysis 
routine computes the inviscid flow past a given airfoil in 
two steps. The region exterior to the airfoil in the 
z-plane is mapped conformally onto the interior of the 
Unit circle in the c-plane, ana the partial differential 
equation (2.9' for (J>(x,y) is expressed in terms of the 
varj..b (r,ai), where ? = r e^‘^. The resulting nonlinear 
equation is then approximated by a finite difference scher.e 

which is solved by a relaxation procedure to provide the 
solution (j> (r,oj) . 

For the inverse problem, we are given the velocity 
distribution Q(s) rather than the coordinates of the 
airfo-1. The basic idea is to use Q(s) to determine both 
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tho mappln9 2 = f(C) of the unit circle onto the desired 
dix^foil and the potential With this free boundary 
approach, the flow calculations are^ done in a fixed 
computational domain and the geometry is determined by 
introducing the appropriate mapping as an additional 
unknown* This results in coupled nonlinear equations for 
<{> and for the mapping function f which we solve iteratively 
using existing routines in the analysis code. 

The iterations go roughly as follows. We start with a 
fj.rst guess for the potential function v/hich we take to be 
the incompressible solution obtained by replacing the 

partial differential equation (2.9) by Laplace's equation. 
..he values of ^ ^ are used in the equation for the mappi ng 
function, which we solve for the approximation 2 = 

Using the mapping in the flow equation then provides 

a better approximation 4>^^^(r,w) to the potential, and 
the process is repeated until the approximations converge. 

In terns of the analysis code, at each cycle the design mode 
starts with an approximation to the desired airfoil 
maps it to the unit circle, and solves for the flow past 
in the usual way. The resulting speed distribution is then 
compared to the desired speed distribution 0(s), and a correc- 
tion to p^*^^ IS made to obtain the now approximation 
The iterations continue until the computed speed distribu- 
tion agrees with the prescribed distribution Q(s) within 
an acceptable accuracy. 
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We now describe this procedure in more detail. Consider 
a conformal mapping z = f(;) from the unit circle onto the 
exterior of an airfoil such as appears in Figure 2. Vie 
assume f has a pole at the origin and takes the point 
C = 1 into the tail of the profile (x(0),y(0)). The 
included angle at the trailing edge will be taken to be 
zero, although the less important case e > 0 could be 
treated similarly. The derivative of the map function has 
an expansion of the form 


(3.1) 


H = f(U 


1 - C r k 

2 “ I ' 

C k =0 ^ 


where the behavior of the mapping at the trailing edge of 
the airfoil is taken into account by the factor (1 - ?) . 

The mapping determines a boundary correspondence 
between the unit circle 4 = and the airfoil which we 
write as s = s(oj), where s is arc length along the profile. 
If we denote the inclination of the tangent to the airfoil 

by 0(s), as in Figure 2 , then on the unit circle c = 
we have 

(3.2) f'(^) = |||-| exp |i arg | 

= M e^Cp {i(0(s(o!)} - 0 ) - 1 
°k ~ ^k ^ (3.1) gives 
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(3.3a) log | = I ‘^cos kw - bj^sin kw 

^2 sin ~ k=0 ^ ^ 


(3.3b) 0(s(o))) + ^ + Ti = I b.cos kw + a. sin kio . 

^ k=0 

Tliese equations show that the mapping is essentially 
determined once the correspondence s = s(u>) is known. 

Under the change of variables z = f(r;)* C = c 
the partial differential equation (2.9) for <}> becomes 


(3.4) r^(c^- - 2ruv?j.^ + (c^- 

+ r(c^- v^)5 + ^ (u^+ v^)[r^uh.+ r^vh ] « 0 , 

2 2 2 2 2 

where s*'(r,w) = 4>(x,y), h = |dz/d^l , u = 4>^/h , 

= (jl^/(r^h^), and c^ is given by Bernouli's law (2.7). 

Note that the mapping function f appears in (3.4) through 
the Jacobian h. Furthermore, for the inverse problem the 
solution Mr,u)) of (3.4) can be used together with the 
data Q(s) to provide boundary values for the determination 
of f'(«;). The relation ^(l,u))= <{> (x [s (m) j ,y (s (w) ) ) yields 
upon differentiation 


(3.5) 


ds _ 1 3^ 

du) Q (s (uiTT *5w 


[l,toJ 


t 


which can be used in the expression (3-3a) for log | f ’ (e^“’) | . 

The problem is therefore described by the equations 
(3.1), (3.3a), (3.5), and (3.4), together with the appropriate 
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boundary cdnditions for to supplement (3.4). The 
iterations used in the computation to solve this problem 
start with the approximation <J>^®^(r,a>) provided by 
incompressible theory. Tf we introduce harmonic function 
G(C) = log (1-C) f‘(C)l# the iteration scheme then 

proceeds by solving in succession 


(3.6) 


ds 


(n) 


d(o 


(to)) 


9 ^ 


(n-1) 


9o) 


(l,w) 


(3.7) 


(C) = 0 


G<”>{e^“) = 


log 


2 sin ^ 


ds 

du) 


(n) 


(m) 


(3.8) 


h^"^ (C) 



exp G^"^ (C) 


t 


(3.9) 


+ ru_v 


n n' "rr 


(n) 
n'n*rtJL) 


+ 


(c2-v2)^ 


(n) 

OJO) 


* r(c^- vliiM . i (uV)Ir\„h^"> . r\h^">) = 0 


for n = 1,2,3,..., where /h*"* , = e V (r^h , 

2 2 7 

and c IS given in terms of u and v by Bernoulli's 
law (2.7). The boundary conditions used to solve (3.9) 
are those of the direct problem. 

The mapping computation (3.7) can be done rapidly using 
the fast Fourier transform to evaluate the coefficients 
appearing in a truncated series expansion for 
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The flow computation is performed using a nonconserva- 
tive difference scheme similar to the one first developed 
by Murman and Cole [27j . Its main feature is type- 
dependent differencing which captures shocks over tv;o mesh 
widths by effectively producing an artificial viscosity 
in the supersonic regions. 

The iterative procedure works very well for subsonic 
flows, presvunably because the initial guess is a good 
appiToxxination to th© solution. In £©ct w© prov© in 
Section 5.2 that a similar iteration converges to a solu- 
tion provided the maximum Mach number in the flow is small 
enough . 

Fo_- the case of transonic flow, the problem is 
complicated by the possible presence of shocks in the flows 
past the various approximations to the desired airfoil. 

A large gradient in the derivative of $^”"^^(l,oj) appearing 
in (3.6) is undesirable, since a discontinuity in (3.3a) 
causes a logarithmic singularity in (3.3b), which is 
inconsistent with the assumed smoothness of the airfoil. 

One way to avoid this difficulty is to solve the equa- 
tions on a coarse mesh. The coefficient of the artificial 
viscosity implicit in the Murman-Cole scheme is of the 
order of a mesh width. If the grid is coarse enough, weak 
shocks are suppressed by this viscosity, as illustrated 
in Figure 4. This smoothing effect allows the process to 
converge even in the case of transonic flow. If the grid 
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is refined, unwanted shocks may appear in the flow, 
causing the iterations to diverge, on the other h 
a solution computea on too coarse a mesh may not accurately 
mode, the actual flow past the airfoil because of the 
smoothing effect of the artificial viscosity. 

in oider to operate on a fine mesh, we have added 
an additional smearing term to inhibit the formation of 
shocks, we describe this term in more detail m the 
aekt section. It has the form of an artificial viscosity 
multiplied by a coefficient s,Aw , where dw is the 
mesh width in the angular direction. The facto, c 
be varied to change the amount of smoothing used. This 
permits the use of more viscosity in the early itera- 
tions when It is important to suppress shocks, and less 
viscosity towards the end of the computation when a more 
accurate solution is desired. The additional smoothing 
term therefore significantly increases the versatility of 
rloeirtn routine* 
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2 . 


The Flow Computation 


In this section we discuss the difference scheme 
used in the flow calculation and also give the form of 
the additional artificial viscosity term used in the 
design procedure. 

For computational purposes it is convenient to 
remove the singularities of ${r,co) and h(r,o)) by defining 

(^(r,o)) = cos (o)+ a - bp) + 4.(r,w) , h(r,w) = 

mv 

The equations for $(r,o)) then become 


(3.10) r (c - u^)<I> - 2ruv<I> h 

rr ro) 


12 2 

+ - (u + v‘) (ruH^ + vHj^] 


(c2-v2)$^^+ r(c2-2u2-v2)4.^ 
0 , 


where 

2 Qq 

u = (r - q^e cos (u+a-bQ) ]/H , 
a- 

V = (r4>^^ - q^e sin(o) + a - bQ)]/H , 

and c2 is given by Bernoulli's law (2.7). The boundary 
conditions (2.13), (2.10), and (2.12) become 


(3.11) 

(3.12) 

(3.13) 


4>(0,(o) 

~ ■ 2? tan(o)+a-bQ)) , 


^0 

= q„e cos (co + a - bp) , 

$^(1,0) 

^0 . , 

- - q„e sin (a - b^) . 
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In the flow computation, centered differences are used 
to approximate the coefficients of , and , 

as well as all of the lower order terms in (3.10). A rotated 
difference scheme due to Jameson [20] is used to evaluate 
the second derivatives . This method uses centered differ- 
ences at all subsonic points. At supersonic points, 
one-sided differences that are retarded in the local stream 
direction are used to produce an artificial viscosity 
similar to (3.15) below. The resulting nonlinear algebraic 
equations are solved using line relaxation in the direction 
of the flow. 

For the two-dimensional flows past a wing section that 
we consider, the direction of supersonic flow is to a good 
approximatxon alligned with the uj— coordinate direction. 

The original Murman-Cole scheme would suggest the approxi- 
mation . 


(3.14) ^ 


It. .-24>. . ,+ 0 - . - 

n 1 , 1-1 1 , 1-2 

(Ao))^ 


which is first order accurate at (i Ar, j Au) . The dominant 
truncation error in (3.14) is Ao)'5^j^^(i Ar, 3 Aw), which has 
the effect of an artificial viscosity. We may consider 
this scheme as an approximation to the equation 

(3.15) r^(c^-u^)4>j.^+ ... = Au) max [0, (v^-c^) ] . 
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The artificial viscosity on the right is absent in the 
subsonic regions and the coefficient tends to zero 
at the sonic line. 

In order to improve the convergence of the design 
routine^ we have added an additional artificial viscosity 
to the flow equation (3.10) . The added term has the form 

(3.16) e.Ao) I— (V(M) (p 1 

which is motivated by the Murman— Cole artificial viscosity 
appearing in (3.15). Here V(M) is a smooth function of 
the local Mach number M which vanishes for M £ and is 
one for M ^ Mj^. We choose the numbers and so that 
this term is effective across the sonic line. For example, 
we may use = 0.85 and = 0.95. This is in contrast 
to the behavior of the artificial viscosity in (3.15), 
which is switched off at the sonic line . 

The term (3.16) is added directly to the rotated 
difference scheme, so that for = 0 the original scheme 
remains unchanged. Figure 5 illustrates the smoothing effect 
of (3.1^) for a flow with a weak shock on a fine mesh. 

Note that the shock does not appear on a cruder mesh . 

By adding the term (3.16) to the partial differential 
equation (3.10) we can obtain satisfactory convergence of 
the design scheme on either crude or fine meshes as desired. 
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3. The Conformal Mapping 


In this section we discuss the calculation of the 
mapping function from (3.3a), and we also explain how the 
Mach number , the coefficient ..f lift , and the 
rotation factor bQ are obtained • 

Durin 9 each design iteration we use the data Q(s) 
and the previous estimate (r,w) for the potential 

function to calculate a boundary correspondence s = s (w) 
between the unit circle and the nth approximation to the 
airfoil. The values (r,o)) are obtained from the 

analysis routine, which uses dimensionless units that are 
normalized by the free stream velocity q„ and the chord 
length L of the profile. Since and L are not specified 
for the design problem, it is necessary to adjust the scal- 
ing at each iteration to make the prescribed data compatible 


with the units used by the analysis routine. 


In order to use the analysis routine we need to supply 
values for the free stream Mach number and the lift 
coefficient C. . To determine we use the relation 

ii 

2 


(3.17) 



Y+1 

Y-1 



1 


which results from Bernoulli's law (2.14) and (2.7). The 

speed q corresponding to the data Q(s) and c* is obtained 

~ ( 0 ) 
Iteratively. In the first cycle we use the value for q,„ 

provided by the incompressible solution. At the nth step, 


34 



• 1 - 

IS chosen to be the scale factor that minimizes the 
expression 


(3.18) 


* - 


,(n) 


where the points are evenly distributed around the unit 
circle and are the velocities along the surface of 

the (n-1) St approximation to the airfoil as computed by 
the analysis routine. 

The coefficient of lift supplied to the analysis 
routine. 


r 

^0 


Q(s')ds 


“ 1 




($] 

1 ~ 

2 


is also affected by the scaling and must be similarly 
adjusted. 

In order to evaluate s^”^(to) we use the data Q{s) and 

( 1 , 0 )) as follows. Formula (2.15) is integrated to 
provide the function 


4> 


l(s) = I Q(s') ds’ . 
0 


h?s a minimum at the stagnation point of Q, say , 

and is monotonic on either side of s„. We define a smooth 
function 


(3.19) 


^2<s) = { 


- - »^(s„) , s < s„ , 

+ ■'•'llaj - «^(S(,1 , e > S 
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which is monotonically increasing ani can be inverted. 

Both the speed Q and the arc length s may then be considered 

functions of the modified potential ^ 2 * 

At the nth stage of the iteration r the potential 
is modified as in (3.19). The result is 
scaled to have the same range as $2 inserted into the 
appropriate expressxons for Q and s. This provides an 
expression s = s^”^a)) which can be differentiated and 
used directly in (3.3a) instead of using (3.5), which 
requires special treatment at the stagnation point. 

The series 


(3.3a) log 


ds 


(n) 


. 0 ) 
sin 2 


do) 


-] = 


I 

k=0 


aj^”^cos 


k(jj 


- 


in ku) 


is truncated at N terms and the coefficients a^ and 

+ ib^"^ , k =1,...,N-1 are obtained using a fast 
k k 

Fourier transform. This procedure does not provide the 
coefficient » which determines the orientation of 

the airfoil with respect to the coordinate axes. To find 
we appeal to the Kutta condition 

0 

(3.13) e sin (a - b^,) . 

At the nth stage of the iteration, bQ ^ is determined so 
that (3.13) will be satisfied as the iterations converge. 

In summary, at each iteration we rescale the data to 
update and , and determine the mapping coefficients 
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a.Q and Cj^ = + ibj, / k = which are used as 

input to the analysis routine. The analysis routine then 
provides the potential ^(r^oi) alon^ with the necessary 
rotation factor to complete the cycle. The iterations 
continue until the velocity along the airfoil computed 
by the analysis code agrees well enough with the prescribed 

j 

data Q(3) . The remaining boundary conditions are automati- 
cally satisfied, since at each iteration we use the analysis 
routine to solve for the flow past a given airfoil. 



4. 


The Boundary Layer Correction 


For the computations to be of practical value it is 
important to supplement the inviscid equations discussed 
thus far with equations describing the flow near the 
surface of the wing section v/here viscous effects cannot 
be disregarded. To do this, the inviscid theory is used 
to design a profile with a finite thickness between the 
upper and lower surfaces at the trailing edge. Next a 
boundary layer correction is computed on the basis of 
the inviscid pressure distribution. The displacement 
thickness of the boundary layer is then subtracted from 
the coordinates of the inviscid profile. Thus the end 
results of the computations are the actual coordinates of 
the airfoil, a viscous boundary layer next to the surface 
of the airfoil, and inviscid flow outside the streamline 
determined by the boundary layer. 

The method oI Nash and Macdonald [28] is used to 
compute the turbulent boundary layer correction. The 
momentum thickness 9 and the displacement thickness 6 
are calculated from the von Karmen momentum equation 


(3.20) 



(2 + H - M^) 


dq Q 
ds q 


T 



where H “ 6*/0* is the shape tactor and t is the skin 
friction. and q are functions of arc length s determined 

by the inviscid solution, and H and t are given by 
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serai -empirical formulas. The ordinary differential equation 
(3.20) is integrated from transition points that must 

be prescribed on the upper and lower surfaces of the airfoil, 
and a starting value for 0* is obtained from the specified 
Reynolds number of the flow. 

Separation of the boundary layer is predicted when the 
Nash-Macdonald parameter 

C = - ~ ^ 
sep “ q ds 

exceeds 0.004. it is important to choose the input speed 
distribution so that remains around 0.003 on the upper 

surface near the trailing edge. This is our version of a 

to Stratford for avoiding boundary layer 
separation [35] . 

When theoretically designed airfoils are evaluated 
in wind tunnel tests, it is sometimes found that the effects 
of the boundary layer cause losses in lift an-' 'ther 
discrepancies between theory and experiment. However, 

with such a Stratford pressure distribu- 
tion using a similar inverse formulation [4] have been 

found to meet their design specifications in wind tunnel 
testing . 
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IV. COMPUTATIONAL RESULTS 


In this chapter we present some results produced 
by the design mode of the analysis code. We include 
a description of pressure distributions that generate 
profiles with no significant drag creep according to 
the analysis code. Possible extensions of the main 
ideas to other problems in transonic flow are discussed 
in Section 4.2. 

1. The Design Procedure 

The inverse method of airfoil design uses as input 
the pressure distribution rather than the airfoil 
coordinates, xn order to obtain airfoils with acceptable 
drag levels, an appropriate pressure distribution must be 
prescribed. In this section we discuss a method of using 
the design mode that produces wing sections with low wave 
drag as predicted by the analysis mode. 

Our first example appeal's in Figures 6 and 7. Figure 3 
shows the speed distribution used to produce the airfoil 
of Figure 6. On the upper surface the speed distribution 
rises from the stagnation point at the nose to a flat 
section of supersonic values along the first sixty percent 
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of chord, and then falls into a Stratford distribution 
near the tail. The distribution over the lower surface 
is entirely subsonic and is arranged so that the Ixft 
is evenly distributed along the section, with aft load- 
ing at the tail. The profile has been provided with a 
gap at trailing edge so that a boundary layer correction 
can bo removed from the displayed coordinates, as shown 
in Figure 7. 

Our experience with the design routine to date 
indicates that drag creep can be avoided by designing the 
to have a small enough supersonic zone. The 
supersonic zone can be increased or decreased in size 
as desired by varying the critical speed c* used in the 
design routine. If too largo a supersonic zone is used 
at design, the middle part of the zone tends to collapse 
at speeds below design, giving rise to one or two shocks 
that can cause significant wave drag at off-design condi- 
tions. This effect is reduced by designing at a lower Mach 
number with a smaller supersonic zone. 

Figure 8 shows an airfoil design with a speed distri- 
bution similar to that of Figure 6 but with the critical 
speed c* lowered so that the supersonic zone is significant- 
ly larger. The pressure distribution was altered slightly 
near the nose and tail of the airfoil to retain the same 
thickness-to-chord ratio and about the same gap at the tail. 
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When drag rise curves are computed for these two airfoils, 
the profile designed with the larger supersonic zone 
exhibits drag creep as illustrated in Figure 9. The 
observed difference in the drag levels for the two profiles 
is due to increases in both the wave drag and the form drag 
of tho second airfoil, although both airfoils have virtually 
identical form drag at subcritical speeds. 

The design mode can also be used to improve the 
performance of airfoils by altering off-design pressure 
distributions. For example, we may exploit the fact that 
some shockless airfoils designed by hodograph procedures 
exhibit characteristic off-design distributions when 
evaluated near the design angle of attack with a lower Mach 
number (cf. [2], p. 96; [3], p. 131, p.l43) . The speed along 
much of the upper surface is roughly sonic, with a pronounced 
peak near the nose of the profile. Such peaky distributions 
also recall the experimental work of Pearcey [ 31] , as well as 
Boerstal and Uijlenhoet [9 ] and Nieuland and Spee [29], 
who have published examples of shockless airfoils designed 
with peaky distributions . 

We illustrate the use of this observation with another 
example. Wo begin with an airfoil that was obtained 
using the design mode with a Mach number = 0.745. We 
use the analysis routine to compute flows past this profile 
with the same angle of attack but with smaller Mach numbers. 
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At = 0.710 there results the distribution shown in 
Figure 10, with an upper surface distribution that 
resembles the characterisiic distributions except for 
a bump in the distribution around sixty percent of chord. 

This distribution is obtained as output from the code in 
the form of punched cards. We modify the distribution by 
removing the bump so that the distribution remains relatively 
flat along this section of the airfoil. The resulting 
distribution is used in the design mode to obtain the 
airfoil shown in Figure 11, in Figure 12 we display drag 
rise curves for this airfoil and a shockless airfoil with 
similar specifications designed by Dr. Jose Sanz using 
the hodograph code of [4] . The airfoil produced by the 
design mode of the analysis code compares quite favorably 
with the shockless airfoil. Figure 13 shows that a near- 
shockless flow is obtained at = 0.740. 

The two previous examples illustrate the procedure we 
use to obtain airfoils with low wave drag. We start with 
an upper surface speed distribution similar to the one 
appearing in Figure 3. This portion of the distribution 
determines the wave drag experienced by the airfoil at 
°^^~^®sign conditions and also determines the growth of 
the boundary layer near the tail. To obtain airfoils with 
a given gap at the tail, thickness-to-chord ratio, and lift, 
the lower surface distribution should be modified as we 
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will indicate in Section 6.1. The value of c^ used 
determines the free stream Mach number M , as well as 
the size of the supersonic zone. In order to find the 
proper size for the supersonic zone, it may be necessary 
to make a tentative choice for c* and use the analysis 
code to calculate a drag rise curve for the resulting 
profile. The size of the supersonic zone should be 
decreased if the wave drag is too high at speeds below 
design. 

The size of the supersonic zone used at design 
determines the height of the pressure peak near the nose 
of the profile at off-design conditions, which in turn is 
related to the amount of wave drag occurring below design. 
The amount of wave drag near the design condition is 
effected by the curvature of the profile at the rear of 
the supersonic zone, which is governed by both the 
prescribed pressure distribution and the amount of 
artificial viscosity used in the design routine. Rather 
than attempting to adjust this area of the profile when 
it is in the most sensitive region of flow, we have found 
it more convenient to make any necessary modifications in 
an additional design run at a lower Mach number correspond- 
ing to the characteristic off-design condition as in 
Figures 8 and 9. To do this, the analysis mode is used 
to obtain the flow past the profile at the design angle of 
attack but with lower Mach numbers . At some Mach number 
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the flow should have a pressure peak near the nose 
followed by a section of nearly constant sonic flow. 

The pressure distribution over the region corresponding 
to the rear of the supersonic zone at the design condition 
IS examined for any irregularities, and, if necessary, the 
pressure distribution is modified near this point, so that 
it more closely resembles the characteristic off-design 
distribution of shockless airfoils. 

In taking this approach, our philosophy is therefore 
to use a relatively simple upper surface distribution as 
in Figure 3 at the design condition. At this stage we 
adjust the remainder of the input distriubtion so that the 
airfoil has the desired specifications and we determine 
the size of the supersonic zone so that the off-design 
performance is acceptable, m doing so we operate on the 
fine mesh of code H with enough added artificial viscosity 
to ensure convergence of the scheme. 

If the .171-Jysis mode is used to evaluate the airfoil 
at the des.^- condition, the resulting pressure distribution 
usually agrees with the assigned pressure distribution 
except near the rear of the supersonic zone where the 
extra artificial viscosity used in the design mode has 
its largest effect. Rather than attempting to achieve 
better agreement betwe-'n cV» gn and analysis in this region 
by using less artificial viscosity in the design mode, we 
instead go to the off-design condition to make any necessary 
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modifications to the profile at this point. Small correc- 
tions usually do not significantly alter the specifications 
of the airfoil that were determined at design. The result 

IS a smooth profile with low wave drag at off-design condi- 
tions . 

We discuss the implementation of this procedure, 
together with the necessary boundary layer correction, 
in Sections 6.1 and 6.2. 

Figure 14 shows an airfoil with a larger supersonic 
zone designed on a fine mesh using a relatively small 
coefficient = 0.05 in the additional artificial 
viscosity term (3.16). The appearance of the sonic line 
suggests the presence of a shock in the interior of the 
flow region which weakens as it approaches the profile. 

This picture illustrates the fact that this approach does 
not produce shockless airfoils, but can provide some contro) 
over the shock strength at the body by fitting a smooth 
pressure distribution. 



2. Extensions of the Technique 

A procedure similar to the one presented here would 
allow the design of transonic cascades with low wave drag. 
Codes whxch compute transonic flow past turbines and 
compressors by using relaxation schemes similar to the 
one in the analysis code used here have been written [19] 
and would presumably lend themselves to a similar design 
modification. An attractive feature of this approach 
would be avoiding the complicated paths of integration 
necessary for the design of transonic cascades using complex 
characteristics in the hodograph plane [4J. For example, 
it might prove possible to obtain cascades with a smaller 
gap-to*“Chord ratio than can be obtained using the hodograph 
method. 

An important extension of this method is to the case 
of three-dimensional transonic flow past wing-body combina- 
tions. The results obtained so far with the design routine 
suggest that by choosing the proper pressure distribution, 
a satisfactory wing might be obtained with the mesh widths 
usually available to three-dimensional codes. Analysis 
codes that compute the transonic flow past a given wing 
are currently available [3/21], It would be necessary 
to extend the method used in two dimensions to treat the 
more complicated free boundary. One possibility would be 
to use a separate conformal mapping at each wing section 
to define the wing. 


47 




V. A CONVERGENCE THEOREM 


This chapter treats some theoretical aspects of the 
design problem. Section 1 describes the simplest design 
problem in incompressible flow, where the method xs exact. 
Section 2 outlines a convergence proof for an iteration 
scheme similar to the one used in the computations. The 
estimates needed for the proof are described in Section 3. 

1- The Incompressible Problem 

In this section the design problem is illustrated for 
the elementary case of incompressible flow. The explicit 
solution obtained here is used in the next section as the 
basis for a convergence proof of an iteration scheme similar 
to the procedure outlined in Chapter 3 in the case of 
subsonic compressible flow. 

To make the presentation as simple as possible we will 
consider the case of purely circulatory flow around a 
smooth object. The Kutta-Jouknwski condition (2.12) is 
then unnecessary and the speed of the flow at infinity is 
zero, which simplifies the asiTnptotic representation (2.13) . 
It is also convenient to formulate the problem in terms of 
the stream function if; as well as the velocity potential ({> . 

We first establish some notation which will be useful 
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in the next section. Consider a flow circulating around 
a smooth body as indicated in Figure 15. We choose units 
so that the total arc length of the body is 2u, so that 
the density tends to one at infinity, and so that the 
maximum speed on the surface of the body is one. We 
measure arc length s from a fixed reference point on the 
body and express the coordinates as functions (x(s),y(s)) 

of s for 0 < s < 2 . The speed of the fluid along the 
body is denoted by 

Q(s) = |u(x(s) ,y(s) ) I 

for 0 < s < 2tr. we let 6 = min Q(s) > 0 so that S < Q(s) < 1. 
The potential function 

s 

fS.2) 4 ,(s) = j g(s') ds' 

0 

is then monotonically increasing and <!> (2 ti) ( 0) = -r > 2 tt 6, 
where r < 0 is the circulation of the flow. 

It is convenient to consider Q(s) as defined by (5.1) 
to be a 2Tr-periodic function defined for -« < s < «», and to 

consider $ to be defined on the whole real line. The 
function 

(5.3) s = S(4>) 

inverse to <Hs) then satisfies 

S(<I> + D = S(<I>) + 2ti 
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for -«o < 4> < » , and the function 


(5.4) 


Q(4>) = q(s(4>)) 


has period -F over the real axis. 

In this section the motion is assumed to be incompres- 
sible, which means we may take p = 1 in (2.2). Formulas 
(2.8a) and (2.8b) then show that the complex function 

(5.3) X(z) = <{'(z) + it|)(z) 


is analytic with derivative 


dx 

^ = u - IV 


Near infinity, the asymptotic form analogous to (2.13) is 


(5.4) 


X(z) 271 ^ 


We normalize x so that 


(5.5) <Mx(s) ,y (s) ) =0. 

We begin by observing that the body is determined up 
to a rotation and translation by the function Q(s) . To see 
this, consider the conformal mapping z = f(c) which takes 
the interior of the unit circle in the ^-plane onto the 
region exterior to the body in the z— plane. VJe assume the 
pole of f is located at C = 0 and f (1) = x(0) + iy(0) . 

In the c-plane, the complex po central x assumes the 
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simple form 



” 2irT 

and in particular for c “ 

x/^iu. Fu) 

He ) = - ^ 

corresponds to the function 4>(s). This provides a corres- 
pondence s = s (tj) betweeri the unit circle and ^he surface 
of the body of the form 

(5.6) = = =(-§). 

which is valid for all w. To determine the body from the 
boundary correspondence (5.6), we note that if 
F(c) = -c^fU) , \ ^ have 


|P(ei")| . I If (e^“) 


ds , , 


and therefore 


(5.7) log|F(e^“)| = log ^ ^ 


d(> do) 


= log 


-r 




This determxnes the boundary values of the harmonic function 
log|F(5)|, If G(c) is a conjugate harmonic function for 
log I F ( C ) I » we have 


(5.8) 


-2 


f ’ (C) = C" |F(c) I exp i{G(U + b } 


where b^ is a real constant. The mapping f(c) is therefore 
determined up to a translation and a rotation by Q(s) . The 
flow 
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(5.9) 


r 


u - iv = 

2TiiC fSC) 

-ib- 
is determined as well up to a multiplicative factor e . 

This is all that can be expected, since a Euclidean 
transformation of the coordinates leaves the speed distri- 
bution Q unchanged. 

We may now change our viewpoint and let the formulas 
(5.7), (5.8), and (5.9) determine a nonzero analytic 
function f'(0 and a flow u - iv(c), say with b^ = 0, 
from a prescribed function Q(s) . Provided that the func- 
tion f(?) determines a reasonable profile, the boundary 
condition (5.7) shows that the resulting flow u-iv(z) 
does have magnitude Q on the body, since 


|u-iv(x(s) ,y(s)) 



<j(x(s) ,y(s)) I 


- 1^1 1^1 
' 'ds' 


= Q(s) 


The question arises whether every smooth, periodic 
function Q(s) determines a reasonable profile (x(s),y(s)j. 
This IS not the case. For example, if 


then 


f • (O 



v!Xp [ C C” 

n=0 " 


f 



body 

where 


I f (C) dc = 2TTi c, e°® , 

| C|=1 
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dui . 


2v 

= i J log|F(0| 


This imposes a compatibility condition 

2n 

0 = I log q(=^Jo-^“ d» 

0 

on Q(s) in order to obtain a closed body in the z-plane: 

In the transonic design problem for flow past an airfoil, 
an analogous condition on Q allows one to design airfoils 
that have a finite thickness between the upper and lower 
surfaces of the trailing edge in order to represent a 
wake extending downstream from the tail of the airfoil. 

A more subtle detail is that the mapping z = f(,-.) 
determined by Q(s) may define a profile with self-intersect- 
ing boundaries. m the airfoil design problem, this 
consideration is important since the body doternuned by 

Q(s) may have so much curvature that the top and bottom 
surfaces overlap. 

Novertholess, we emphasize that the formulas (5.7), 
(5.8), and (5.9) do provide a locally one-one mapping 
2 = f(0 and a flow u-iv(r.) if only 0(s) is positive and 
periodic. Similarly, the numerical computations for the 
transonic design problem are possible under very mild 
requirements on Q(s), and the process converges whether 
or not the resulting airfoil is physically realizable. 


53 



It falls to the user of our computer code to make the 
modifications of Q(s) necessary to obtain an acceptable 
geometry • 
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2. 


We now consider the more compUcated case of subsonic 
compressible flow aroend a smooth obstacle. We wish to 
Show that an iterative procedure similar to the one used 
in the numerical computation converges to a solution of 
the inverse problem. The approach we take exploits the 
idea that incompressible flow can be considered to be 
e limiting case of compressible flow as the speed of 
sound c becomes infinitely large. The convergence proof 
requires the maximum Mach number in the flow to be small, 
which can be assured by taking the prescribed critical 
speed c. large enough. In this case we obtain a Poisson 
problem with nonlinear inhomogeneous terms which we solve 
by Iteration. We are able to use standard estimates 
e.xpressed in terms of Holder continuity to show flat the 
iterative procedure defines a contraction, so that the 
iterations converge. Less restrictive results could prob- 
ably be obtained using deeper techniques from the mathemati- 
cal theory of subsonic flows. However, the proof outlined 
here is a satisfactory illustration of the computational 
procedure, which is our main concern. 

For the subsonic design problem, the critical speed 
c* is given in addition to the speed distribution Q(s). 

We continue to use the conventions of Section 5.1, with 
units chosen so that p approaches one at infinity and 
max Q(s) = 1. As in the computational procedure, we solve 
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the free boundary problem by finding both the map z = f(C) 
from the interior of the unit circle to the region exterior 
to the desired body and by solving for the compressible 
flow u(C). We show that for a fixed speed distribution 
Q(s), if c* is sufficiently large we may construct a map 
2 = f(C) and a streeun function both depending on c* , 

which approach the corresponding incompressible solutions 
determined by Q(s) as c* -»■ » . 

For compressible flow, the velocity potential (|> and 
stream function i(», considered as functions of the variable 
C = C + in» satisfy 


(5.10a) P<P^ — , 

(5.10b) P(j)^ = , 


We choose to work with t{) instead of ({) so that we may use 
the Dirichiet boundary condition i/) = 0 rather than the 
Neumann condition d<p/dv = 0. 

By eliminating <}. from (5.10a) and (5.10b) we obtain 
the equation 


(5.11) ^ 

c 


^ 


where u^ = «{»^Vp^|f'|^ , v^ = u^+ v^, 

2 2 Y-1 

and c = (Y+1 )c^P /2. The density o can be obtained from 

Bernoulli's law in the form 
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(5.12) 






2p^|f (C) 1 


l(Y-^l) 

2TFIT 


2.Y-1 


c*P 


1 Y+1 

2 Y=r 


TliG solu'ti.on we shall obtain has the asyir»ptotic fonu 

'f' -§^r logUI 

as C -»• 0, where the circulation r is given by 

27T 

r = - I Q(s) ds . 

0 

We have 2ir > -r ^ 2tt 6, with 6 = min Q(s) > 0. 

We again consider the analytic function F(c) = 
which satisfies the by now familiar boundary condition 


(5.13) log|F(e^“)| = log 


Q(<{>(e^'^)) 


9_ 

do) 


d>(e^“) 


Using (5.10) we have p 3({)/9w = -diji/dr , so that we may express 
the potential functxon on the boundary in the form 

0) 

(5.14) 4>{w) = - f (e^^^) da, . 

If we knew the solution i|/{re^‘^), (5.14) could be used in 
the boundary condition (5.13) for the determxnation of |f'(c)I. 
We consider the expressions 


(5.15) V(c) = iHc) - (r/2ir) log|r,I , 

(5.16) H(c) = log|F(c)| - log|FQ(c)I, 
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where 


-C derivative of the conformal mapping 

obtained by solving the incompressible problem with the 
same data Q(s) , as described in Section 5.1. y and H are 
perturbations of the corresponding incompressible solutions. 
Note that T and H are not singular at the point C = 0, in 
contrast to both the mapping f(c) and the stream function 
4>(C). 

Substitution of the expressions (5.13) and (5.14) into 
the equations (5.11) and (5.13) for il> and F shows that 
4' and H must satisfy equations of the form 


(5.17) 


(5.18) 


AY = M[Y,H,c*] 

^( 6 ^*^) = 0 
’ AH = 0 

. H(e^“) = N[Y,H,c*l , 


where M and N depend nonlinearly on Y and H and formally 
2 

tend to zero as c* -► •«. We give the explicit form of M 
and N in the next section. M is a function of Y and the 
partial derivatives of Y up to second order, H and log|rQ| 
and their first order derivatives, and the variables 
C and n . N involves the boundary values of H,V , and 
the normal derivative (SY/Or) (e^^"') in a manner 
similar to (5.13) and (5.14). Tjie function Q(s) enters 
the problem through the term N[Y,H,c*}. 
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We are interested in solving equations (5.16) and 
(5.17) by iteration. We set = 0, and formally 

define an,^ ||(n+l) solutions of 

(5.19) I 

'p(n+l)(ei'^) = 0 
’ = 0 

(5.20) 

.H(n+l)(ei“) = N . 

We denote this operation by 

(5.21) (y(n«)_„(ntl)) = i,(, (n) _„(n) _ 

We wish to show that for a given Q(s), if c* is large 
enough the iterates satisfy an inequality of the form 

(5.22) J-i.^«+l).y(n)B ,H(n+l)_H(n)B < g (n)_^ (n-1) , 

+ J} 

with 0 < 1, which implies convergence of the iteration. 

We therefore must examine how the solutions and 

of (5.19) and (5.20) depend on the functions y (*^^ and 
and we must define the norms to be used in (5.22) . it 
turns out that since M and N are formally small as c* 
we may succeed by using basic estimates for Laplace's 
equation which are expressed in terms of Hfllder continuity. 
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Let D denote the closed unit disk and consider the 
space consisting of n-times continuously differ- 

entiable functions u defined in D with finite norm 


i, j<n 


+ sup 
i+j =n 


u(C2>l 


ICi - C2I 


where n is a nonnegative integer and 0 < a < 1. 
is a complete space with this norm. In addition, if 
u e then u € when n* <_ n 

andf<a, and i '“'n+a- “l '“2 ® <=n+a *“> ' 

than U^-Uj s V„(D) and < K„+„« di'n+a'ttj' n+a ' 

If u e C^^^(D) and G is (n+ 1 ) -times continuously differ- 
entiable on the real axis, then g(u(^)) e C^^^(D). 

We also need the idea of Hdlder continuous boundary 
values. If g(w) is an n-times continuously differentiable, 
2TT-periodic function defined for all 10 , we v/ill say 


g 6 norm 


Igl, 


U) 

i<n 


n+a, 3D 
is finite. The norm !•! 


= sup 1 3^g ((d) I -f sup 


!*"9((Di)-3;;g(,D2)| 


n+OL, 8 D 


a)j^,o)2 - 0 ) 2 !“ 


has properties similar to 


"‘“n+a* ^ f(e^“) e 

and llf(e^^)l . ^rv < K' IfO . 

' n+a,9D — n+a n+a 
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For any fixed a, 0 < a < 1, we will show y 6 q /ni 

(n) ' 

and H e by using the following fact, which is 

an elementary instance of the more general a priori esti- 

estimates of Schauder [ 6 ] ; 

THEOREM. Let u(c) be the solution of 

Au = f 
u(e^“) = .{, 

where f e c^(d) and <p e Cj^^Od). Then u e C 2 +„(D) and 

(5.23) '“*o+„ < ) . 

2+a - S'- a ^ 2+a,dD^ ' 

where depends only on a. 

We will also use the following consequence of (5.23). 

COROLLARY. Let u(?) be the solution of 

’ Au = 0 
, G(e^“) = 4 , 

where ^ e Cj^_j^^(3d). Then u e (D) and 

(5.24) •ul,. < K.lSl 

1+a - ^^4 l+a,9D ' 

where depends only on a. 

This result can be obtained from the previous th-crcn 
(5.23) by setting 

= Izr +(§7) j MW) do)' , 

0 
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where u is obtained by solving 


’ Au = 0 

u 2n 

u(e^“) = J Hu’) du’ - ^ I Hu') du* . 
0 0 


We assume the data Q(s) is three times continuously 
differentiable. The expressions (5.2), (5.3), and (5.4) 
show that the function Q(<J) is also that smooth, and 
we set 


(5.25) 


^ — co<$<oo 


sup 

o<$<o 

3 = 1 , 2, 3 


Q(<t) 


Using (5.24) we see that the harmonic function log|FQ(c)| 
determined by (5.7) is in (D) with I log (F q (-) 1 1 
1 , where Kj d. -.pends only on , 6 , and a. This 

implies |FQ(i;){' (D) with a similar bound 

(r^^,r2) 


0"” ‘ ‘l-^a - ^5 

Consider the closed subset B 


defined as 




2+tx - ^ 1 ' - ^2 


} 


In order to show the iteration scheme (5.21) converges, 
we first establish the following 

LEMMA A. There is a constant depending only on 
a, Y, 6 = min Q(s), and such that for c* > , the 
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operator L(-,.,c*) formally defined by (5.19) and (5.20) 

Is well defined on end maps into itself. 

Thus the iterates and h'"^D all exist and satisfy 

fta 'ita 11. 

We sketch the derivation of the estimates necessary for 
the proof of Lemma A in the next section. There we show that if 

I'^l^ ®(6/2,l) large enough, we have 

Mfy^sH^scJ € c^(D) and Nt'i'^,H^,c*] e C^^^Od) , with 
- K,/c, 2 and 

- ' where Kg depends only on 

a,y,S, and Kq. Using the basic estimates (5.23) and (5.24), 

we tnen have that the functions = LCF^, 

defined by solving (5.15) and (5.16) satisfy 


(5.26) 

(5.27) 


*"'2''2+a 1 


< K /c 
z 1+a — 


with depending on a,y,6, and , which establis..es 
Lemma A with = K_/5. 

A 7 

Lemma A shows that the functions y and can 

indeed be generated for n = 1,2,... . m order to show 
convergence, we have 


LEMMA B. There is a constant depending only on 
a,Y,5, and such that for c* > Kg , the operator 
L(-,.,cJ defined on B^^/^,!) “ a contraction. More 
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specifically, if ~ ^ 

(’i'^jH^) = we have 

‘'*' 2 " ’^4‘2+a'*' '” 2 “ ” 4 *1+0 - ®(*'‘'l‘''‘'3*2+o'^*”l"®3*l+o^ 

where 0 < 1. Thus the iterates 4 '^"^ and form Cauchy 

sequences in the complete spaces and 

To establish Lemma B we consider the expressions 


(5.28) 


A('? 2 -'!' 4 ) = ,c*] - Mr 4 ' 3 ,H 3 ,c*l 

’l'2(e^“) - '{'4(3^^^) = 0 


(5.29) 


A(H2~H^) = 0 


H 2 (e^“) - H^(e^“) = N [4'2/Hj^,c*] - N[Y^,H 3 ,c*] 


In the next section we see that 




3'“3'''*' a 


1 Kg(iy2-y^i2_j,^+ ) , 

where Kg dpends on a,y, 6 , and Kq. Applying the basic 
estimates (5.23) and (5.24) to equations (5.28) 
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and (5.29) 



we have 


'V''4'2« i (■'lO''=*^)('V-3'2l-„+ '«l-«3'i+„) 




whxch together yield the statement of Lemma B with 
" ^10 ® 

The Cauchy sequences of iterates and there- 

fore converge to functions y 6 C 2 ^^(D) and H € 
which must satisfy (-F^h) = L(f,H,c*) if c* > K^, We 
note rhat since L is a contraction, and H are the 
only solutions of (5.16) and (5.17) with J'l'l, <6^2 

and < i. Moreover, the expressions (5.26) and (5.27) 

show that the solution satisfies 

»'J'«2+cc + < 2K^/c*2 . 

Recalling that f and H are perturbation quantities represent- 
ing the difference between the compressible and incompressible 
solutions for a given Q(s), we see that the compressible 
solution indeed tends to the corresponding incompressible 
solution as the critical speed c* - This fact, together 
with the explicit solution available in the incompressible 

case, is the underlying basis of .:he above convergence 
proof. 
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Inequalities for the Convergence Theorem 


In this section we consider in more detail the 
inhomogeneous terms M['?,H,c*] and N['?,Hj^,c^] appearing 
in equations (5.17) and (5.18) and indicate how the 
estimates mentioned in the discussions of Lemmas A and B 
are derived. 

A calculation shows that the inhomogeneous term 
M['i',H,c*] of (5.17) has the form 


(5.28) M[’l',.I,c*l 




_2 Y+l|p I 2 

C*P KqI 


,H, acglFgU ) 
^i ” ^i 

exp 2H 


where = C/ = n» and T is a third degree polynomial 
in its arguments with coefficients depending only on the 
circulation P. The term N['l',H,c*], takes the form 


(5.29) Nl4',H,c*l = log[l + ^ (e^“)J 

- log p (e^“)-log IQ(<I' (u) )/Q(-ru)/2Ti) 1 

Here the density p can be defined by Bernoulli's law (5.12) 
to be a function p = R(q /c*) , where q = (t|;^ + i{'j^)/lf'|“ , 
valid for C 1 1 2Kj^c* , with R(0) =1. R is the 

subsonic branch of the multiple-valued function giving 
p in terras of the gradient of the stream function. The 
constant depends on y, 2K^^ = [2/ (y+ 1) ] . The 

only information about R that we need is that for 
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~2 2 

*3 1 t there is a constant K 2 depending on y such 

that 1 > R > K 2 “^ > 0 and |R'MR"| < K 2 . Recalling the 
expressions (5.15) and (5.16) giving ’i' and H in terms 
of 41 and f • , we have that 


(5.30) 


1 

iFol^exp 2H 



(C'i'^+n4’^) 




The term $(a>) appearing in the expression (5.28) 
for N[4',H,c*J is written in the form 


0) 


(5.31) 4>(u) = (r/c) f i 

J n 


P(e^“) 


_r ^ 3’?(e^“) 

27t 3r 


doj 


rather than (5.14), where the constant c is given by 
2 ti 


(5.32) 


= - J ~ 

J ntn 


P(e^“) 


__r . 3’l'(e^^) 

*< A 


2tt ^ ar 


du) . 


The factor c is introduced to make sure h> (a)+ 27 r ) -0 (o))= -r, 

/N 

SO that Q( 4 >(o))) is 2 Tr-periodic and smooth. 

We first describe the inequalities used in the proof 
of Lemma A. Since the expressions (5.28), (5.29), (5.30), 
and (5.31) are rather complicated, we do not attempt to 
present all the details. Vie have chosen the data Q(s) to 
be smooth enough so that nothing more sophisticated than 
the mean value theorem is needed. 
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we eseume «■,« e . „e claim that if c. is 

large enough, K[g,H,c*] € C (D) with IM['P,H,c*JI < K /c ^ 

and N[y,H,c*] e with I . . < 

2 l+ot^dD — 

^ In the following expressions, K. 

will denote constants depending only on a,6,y, and K^. 

Consideration of the expression (5.30) for shows 
that g^ € with < K^. Recalling the 

properties of the function p = R(qVcJ) defined by (5.12), 
we see that p € c^^^(d) with , provided 

that we have chosen cj > k^'/K^. This choice keeps qVc*^ 
on the subsonic branch of the density-stream function 
relation. We then see from (5.28) that M['i',H,c*] e c ^D) 
and we obtain t.M(’i',H,c*] I ^ . 

We next observe that each of the three terms on the 
right hand side of the expression (5.29) for N[4',H,c*] are 

first term is well behaved since we have 
|r| > 2ir5 and •'*'• 2+0 1 ^/2, giving | (2Tr/D (3 V/3r) (e^“) | _< i 
and its norm is bounded by a constant times *'{'• 2 + • The 
second term log p(e^“) is also well behaved since p is 
bounded away from zero, and we can estimate 


• log P(e 

The desired factor 1/c*^ is essentially due to the fact that 


log p(e^“) = log R(q^ (e^“) (c*^) _ i^g r(q) 


= a!(£^ 


^ R'(tqVc^2) 


I 


R(tqVc*2) 


dt 
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by the mean value theorem. 

The last term in (5.2a) can be estimated Ly 

( 5 . 31 ) Hog Q(0(a,))-iog 0 (-r„/ 2 ,) P , 1/^2 ^ 

This inequality is more complicated and we sketch it as 
follows, using the fact that 0 is three - Imes continuously 

differentiable one car easily obtain 
Hog <5(4>(oj)) - log Q(~ru/ 27 tll, . 

■' l+a, 8 D 

We then express c(w) + rw/2r in the form 

0> 

^(w) + ru/2n = -Il£ 


/ p [l^r + <3u) 

0 

-I-^[lr^llJo»-/||dw 


° d 

and estimate the three expressions on the right in terms 

of r+o, p-l, and y, respectively. The factor r+c can be 
written 


r+c = 


r t r 

= - I fJL + ill f d'P 

We again gam a factor l/c.^ from the expression p-l, and 
we may obtain 

* '■“''^"'uu.sD i KisU/o.^ * iyp„J , 
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and (5.31) follows. We therefore have an inequality 

2 

as stated in Lemma A, provided ^ K^/Kj^. 

We now consider Lemma B. We wish to verify that if 

(4'i,Hi), (4*2»H2) ^ ®(6/2 1^ large enough, we 

have 


(5.32) - M('i'2,H2,c*]l^ 

1 (Kg/c* ) 2+a *^l“^3*l+a^ ' 

(5.33) - N[H-2,H2,c*]I^^^^ 


Consider the first inequality. The form of (5,28) shows 
that, with an obvious notation, we may write 


■«1- «2'« i C''i-’'2'2+«* '"i-«2'u<. 'Vr - ^r' 


+ lexp(-2Hj^) - oxp(-2H2)>^) • 


It is easily seen that the last term is dominated by some 

constant times ^2^1+a* Kxanination of the expression 

-2 

(5,30) for q shows that v/e also have 
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'I , < 

1+a ~ 


K 


2o('V’2'2.<. 



and using this result with p = R(5^/c,2) then gives the 
first inequality (5.32). 

To obtain the estimate (5.33) we wrxte 


*^l“^2“l+rt,DD - ‘log(l+{^Vr)3H'j^/3r (e^“) ) 

- log(l+(27r/D Dy /3r (e^*^))l, 

^ ' 1+a, 3D 

+llog P,(e-)- 

+llog 6(<r_^(io)) - log Q(':2(‘^))‘i+et,3D 

and estimate each of the terms on the right separately. 

The last expression is the most complicated, and can be 
treatea in the same fashion as was the term 
'log Q(<'(u.)) - log in Lemma A. 

Straightforward calculations then show that an estimate 
of the form (5.33) holds, and the statements of Lemma B 
follow. The map (i^^H^) = L (T^ , , c* ) is a contraction 
and has a unique fixed point that can bo obtained by 
Iteration. 

There remains a technical point duo to the introduction 
of the factor c xn the expression (5.30) for 4 ( 0 ). 

We desire the fixed point (v,ll) of the mapping to provide 
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solutions iI'CO and F(C) = f’(C) to the equations (5.11) 

and (5.13), where (5.14) rather than (5.30) is the expres- 
sion for 4 (w) . Therefore we should check that for the 
fixed point (Y,!!) we have 


c 


I 


1 fL_ + LI 

p | 2ti 3r 


dio 


r 


To see this, note that the corresponding function 
i|'(C) = (r/2Ti) logUI + 'P(C) by construction satisfies 


the equation 


+ (4'r,/p)r, ' 


with 4»(c) r/2TT log |^| as ^ 0, whether or not c = - F. 

By Green's theorem we therefore have 


2-n 

- 1- 


27T 


P(e ) 


3 i{^ 
3r 


(e ) 


du) = Ixm 


P (ce ) 


||(cei“)c 


do3 


r 

" " p(0) * 

The expression (5 •30) for q shows that 

p(0) = R(q^(0)/c*^) = 1, so c = - r as desired. This 

completes the convergence proof for the subs->nic inverse 

problem. 
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VI. DESCRIPTION OF THE CODE 


In this chapter we explain how to modify the input 
speed distribution in order to obtain airfoils with given 
specifications. We also describe the other input para- 
meters necessary for the operation of the design mode. 

1 • Achieving Design Specifications 

We refer c gain to the speed distribution illustrated 

. . 3 

in Figure 2. The form of the upper surface distribution 
determines the amount of wave drag experienced by the 
airfoil and the growth of the boundary layer along the 
upper surface. We suggest using a distribution with 
constant supersonic values over the first portion of the 
profile, followed by decreasing values along the rest 
of the surface in accordance with the Stratford criterion 
^seo 0-003 near the trailing edge. 

The value of c* should be determined so that the 
supersonic zone has the proper size, as discussed in 
Section 4.1. 

The lift of the airfoil is related to the area between 
the upper and lower surface speed distributions. By varying 
the lower surface distribution, the lift can be distributed 
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as desired over the airfoil. The free stream Mach number 
is also affected by changes in the speed distribution. 

For example, increasing the magnitude of the velocities 
along the lower surface will generally decrease the lift 

and increase M . 

00 

The thickness-to-chord ratio of the wing section can 
be adjusted by varying the slope of the speed distribution 
near the stagnation point Q = 0. Increasing the slope 
will result in a thinner profile with little change in 
the lift of the airfoil. The free stream Mach number also 
increases as the thickness-to-chord ratio is decreased. 

The vertical separation between the uppi-r and lower surfaces 
is decreased as the slope is increased . 

The relative position of the upper and lov/er surfaces 
at the trailing edge can be adjusted by changing the 
velocity distribution near the tail. The vertical separa- 
tion is Increased by raising the prescribed speed at the 
tail on both the upper and lower surfaces. This wi"*! not 
have a strong effect on the thickness-to-chord ratio. For 
most purposes the vertical separation should be around 
0.015 so that after removing the boundary layer a gap of 
around 0.007 remains. The horizontal separation can be 
adjusted by changing the amount of arc length near the 
tail. 

Finally, we mention that a decrease in the prescribed 
critical speed c* will generally increase the size of the 
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supersonic zone, increase the free stream Mach number, 
decrease the thickness-to-chord ratio, and increase 
the vertical separation at the trailing edge. It should 
have ’j.ttlo effect on rhe lift or horizontal separation 
at the tail. 

When designing an airfoil with given specifications 
it is advisable to proceed in stages by modifying the 
pressure distribution in the appropriate areas one at a 
time so that the effects of each change can be isolated. 
When beginning a new design it is useful to make the 
initial runs on a coarse mesh of 40>:7 or 80x5 5 points, 
where results can be obtained in one or two minutes 
on the CDC 6600. The finer mesh can then be used to 
make minor adjustments to the pressure distribution and 
obtain accurate resolution for the final runs. 
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2. Operation of The Code 

We have included the design modification to program H 
in such a way that when design parameters are not explicitly 
specified, they assume default values that do not affect 
the operation of the analysis routine. The description 
of the operation of the analysis code given in 14] therefore 
remains in effect for the nev7 version also. We assume 

the user is familiar with this description of the analysis 
routine. 

For operation of the design mode, the input speed 
distribution should be provided on TAPE6 in the format 
given in Table 7.1, with negative values along the lower 
surface followed by positive values on the upper surface, 

as in Figure 3. The arc length s must be monotonically 
increasing . 

Table 7.2 contains the other input parameters necessary 
for operation in the design mode. These parameters have 
default values as indicated and can also be specified using 
standard namelist conventions. We provide some sample 
commands below to illustrate their use. 

The parameter NDES specifies the number of overall 
iterations to be performed in the design mode. Each itera- 
tion consists of NS cycles of flow computations, followed 
by a new mapping to the unit circle as described in Sec- 
tion 3.3. Since the time required to calculate a new 


i 
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mapping function is small compared to the time required 
the flow past a profile, we use relatively few 
cycles of flow computation between each mapping. We 
generally specify NS = 10, NFAST = 0, and NREXLAX = 1, 
so that 10 relaxation sv;eeps are used betv/een each 
mapping. The results we have prese.ited have all been 
produced with this choice. 

The parameter T5TEP is a relaxation factor for the 
Fourier coefficients determining the mapping function. 

The rate of convergence of the overall iteration procedure 
depends on the value of TSTEP. if TSTEP is too large, the 
coordinates of the profile may oscillate and the scheme 
converge slowly or not at all. When a small amount 
of artificial viscosity is being used in the flow equations, 
it is sometimes necessary to use smaller values of TSTEP 
to avoid abrupt changes in the successive profiles that 
may cause undesirable shock formation. We have TSTEP 'v 0,2 
^cr the design procedure outlined in Section 4 . 1 . 

An example of the control cards and data cards used 
to design an airfoil on the CIMS CDC 6600 is given below. 

To improve turnaround time we have found it convenient to 
break up the 30 b into several shorter jobs rather than one 
longer 30 b. 

The relevant control cards have the form 
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GETPF (TAPE6=SPEED) 


Speed contains Q(s) as in Table 7.1 


GETPF (LGO=HDES) 


LGO. 

SAVE (TAPE3-COMP1) 


SAVE (TAPE4=COORDl) 


HDES IS a compiled version of the 
modified }:analysis code H. 

Tapes contains data to continue 
the computation if desired. 

Tape 4 contains the coordinates 
of the resulting airfoil in FSYM=1.0 
format, and the final pressure 
distribution. 


For an initial run on a crude grid, the first card 
used as input to the program could be 

( SP RN = 0., ALP=0., NDES=1 $ ] 

NDES - 1 puts the code in the design mode. The specified 
angle of attack with respect to the x and y axis must also 
be given on this card. The Reynolds number is zero for the 
inviscid flow calculations. The speed distribution is read 
from TAPE6 and a plot of Q(s) is provided. (if CSTAR < 0, 
the program will then terminate.) 

The next cards are 

I ?P NS=-1, ITYP=1 $ ] 

The grid is coarsened from 160x30 to 80 x 15 . 

( $P NS=-1, ITYP=1 $ ] 
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The gird is coarsened from 30x75 to 40^*7. 


[ $P NDES=20, NS=10, NFAST=0, 

NRELAX=1, EPS1=0.5, TSTEP=0.2, 

REM=0.5, ITYP=4, XP=1.0, 

KDES=10 ] 

Twi—ity design cycles are performed. After er.ch increment 
of 10 (KDES) cycles, the coordinates of the resulting air- 
foil, the Mach number diagram, and a Calcomp plot of 
the flow are produced. XP=1. causes the desiredyi pressure 
distribution to be plotted on the graph also; if XP=0. 
it will not appear. The results of this computation are 
automatically saved on TAPE3 after NDES cycles . 

[ $P NS=1, ITYP=-1 $ ] 

The grid is refined from 40x7 to 80x15. 

[ $P NDES=20, NS=10, ITYP=4 ] 

Twenty more design runs are performed on the new grid, and 
the results rev;ritten on TAPE3. 

[ $P ITYP=0, XP=0. ] 

The computation stops. XP=0- causes the airfoil coordinate 
to be written on TAPE4, This run takes about 90 CP sec- 
onds execution time. 

If the airfoil produced by this run does not meet the 
desired design specifications, Q(s) is changed and the run 



is repeated. If the airfoil is satisfactory, the run may be 
continued on a finer mesh as follows. 

GETPF(TAPE3 = COMPl) 

GETPF(LOG = HDES) 

LGO. 

SAVE(TAPE3 = C0MP2) 

SAVE(TAPE4 = C00RD2) 

The data cards are then 

{ $P NDES = 1, RN = 0., ALP = 0. $] 

( $P NS = -1, ITYP = 1 $ ) 

Coarsen mesh from 160 30 to 80 15. 

f $P NS = 0, ITYP = 1 $ J 

Read in data stored on TAPE3 (stored v/ith Mxn = 80x15) 
to continue the computation. 

[ $P NS=1, ITYP= -1 $ ] 

( $P NDES - 10, NS = 10, ESPl =1., ITYP = 4 $ ] 

[ $P ITYP = 0, XP = 0. $ ] 

This run takes about 130 CP seconds execution time. 

The next run might have 

( $P NDES = 1 , RN = 0 . , ALP = 0 . $ J 

[ $P NS = 0, ITYP = 1 ] 

[ ?P NDES = 10, NS = 10, EPSl = 0.75, ITYP = 4 $ ] 

[ SP ITYP = 0, XP = 0. $ ] , 

and so on. 
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If It is desired to do the design in a single run, 
the data cards might read: 

I $P NDES = 1, RN = 0., ALP = 0. ^ ] 

{ $P NS = -1, ITYP = 1 ] 

[ $P NS = -1, ITYP = 1 ] 

( $P NDES = 20, NS = 10, NFAST = 0, 

NRELAX = 1, CPSl =0.5, TSTEP =0.2, 

RI^M = 0.5, ITYP = 4, KDES = 10, XP = 1. $ ] 

I $P NS = 1, ITYP = -1 § ] 
f $P NDES = 20, NS = 10, ITYP = 4 $ J 
[ $P NS = 1, ITYP = -1 $ ] 

[ $P NDES = 10, NS = 10, EPSl = 1.0, ITYP = 4 $ J 

t $P NDES = 10, NS = 10, EPSl = 0.75, ITYP = 4 $ J 

( $P NDES = 10, NS = 10, EPSl = 0.50, ITYP = 4 $ J 

( $P ITYP =0, XP = 0. J 

This run takes about 520 CP seconds execution tine on the 
CIMS CDC 6600. 

For the present version of the code, a separate run 
IS required to perform a boundary layer correction for 
the airfoil designed with invlscid the ry. Assuming the 
coordinates and pressure distribution from the design run 
wore stored on TAPE4 = COORDl, the control cards for a 
boundary layer correction would be 
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GETPF(TAPE3 = COORD J ) 

GETPF(LGO = HDES) 

LGO. 

SAVE (TAPE 3 = COORD 2) 

The required data -ards ha^^e the form 

[ $P RN = 20.S06, XP = - 1 ., PCH = 0.07, PLTSZ=8 0 $ ] 

C $P ITYP = 0, XP = 0. ] 

COORD2 then contains the corrected coordinates in FSYM = 2.0 

format. A plot of the profile is also generated by this 
run. 

Finally, we mention that the parametei XOUT can be 
used to obtain speed distributions for use in the design 
mode . The command 

I SP XOUT = 1.0 ] 

will cause the speed distribution currently in memory to be 
written on TAPEJ in the format sliown in Table 7.1. 
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TABLE 7.1. 


TAPE6 INPUT SPEED DISTRIBUTION Q(s) 


^''\Cols. 

CTards 

1-10 

11-20 

P 21-40 

1 

XIN 

CSTAR 


2 

-joj at tail 

initial value of 
arclength 

3 

speed along profile 

1 increasing values 
1 of arclength 

♦ 

• 

• 



XIN + 1 

1q| at tail 

final value of 
arclength 


88 



TABLE 7.2. DESIGN PARAMETERS 
Glossary of Input Parameters for Design Mode 


Parameter Default 


CSTAR 


KDES 


EPSl 


NDES 


PLTSZ 


QPL 


QPU 


REM 


TSTEP 


XOUT 


XIN 


100 . 


10 


- 1 


50. 


.85 


.95 


0 . 


.2 


0 . 


None 


Description 

Real. Critical speed c*. If c* < 0, the 
program plots the prescribed speed distri- 
bution and halts. 

Integer. Graphs and flow printout are 
generated every KDES design cycles. 

Real. Artificial viscosity coefficient 
appearing in the expression (3.16). 

Integer. Number of design iterations 
to be performed. 

Length in inches of profile in 
graph generated when boundary layer 
correction is performed. 

Real. Lower limit Mq of the cutoff 
function V(M) in formula (3.16) 

Real. Upper limit M^ of the cutoff 
function V(M) in formula (3.16) 

Real. Relaxation parameter for determina- 
tion of M . 

oo 

Real. Relaxation parameter for coeffici- 
ents of mapping function. 

Rssl* XOUT 1 causes the current 
'^^locity distribution to be written 
on TAPE3 in the format of Table 7.1. 

The computation is then terminated 

Real. Number of points used in prescribing 

input speed distribution Q(s). 
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LISTING UF COUE 


lQUTPUT!!rA|lE&'*!J^NP0?'?AJto) ' » tOO,TAPE^ - ^Gu, UPE2 ■ 

CJM,1CN/FL/FLUXT^»CD4»COW» IM-CL 

COMMON Prii(l&2>31 ),FP(io2^ 3i ), AOIJ, 6 (31)*C{31),D(31J,E(31) 

CLmMMUN /a/ pi #TP^ fif AD^ E'i^ ALP^F N/ PCHj XPj Tc.run. r pm t n Drt v 
1 #XA# YA» TL#OT»OR#DfcLTri»DcL(<#FA# DCN# DSN# RAA» E PS I L# OCRI T#C I# C2 

I #NPrS#L!lV!pp 2 a K^^''-’’^*'^‘^'^‘^^”''^2,N3,N4#NT,iAX 

5 >SCAL0a#N6#GANPA#NQPT#CSTAR#RfcH#0bP#QINF#TSTEP.XDLT 

6 >INC#QFAC#6a.'1#KDL5#FLTS2#GPL#0PU ^ ^ ^^P# XULT 
DIMENSION C0.')C(67)#CLA(2)#NAMEFR{6) 

EOUl VALENCE ( COhC ( 1 ), PI ),( CL X# C LA ; 1 )),( AL PX#CLA ( 2 )) 
EXTEPNAl\sT£RR^‘^'^'^^”'^^ PRCCESS a NAIELIST ERkiR 
♦**NtlN-ANSI**<= 

^^^*®^^^^®^P*CL#LH#FSYM»GAM.'1A»IS#ITYPtl2«KP»ii i<»pp 

1 ,56/6/^^ ' NANEkP/t,*0/ , 01,02#SL/3*C./ # OPl/.W #XPF/J1J 

AA(1) >. 99999. 

IN0CC=0 
NEW = 1 
NFAST-1 
NRELAX=6 

CALL SYSriNC( 66#NAP.Ekk) 

MA = NA 
REWIND NA 
WRITE (N2»180) 
k£AD (N6»P) 

IF (CL.NE.IOO.) MODL > 0 
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IF (1Z.GE.80) = N2 

IF (NS.Eti.O) Gu TQ 30 
CALL RESTPT 
ClX « LL 
ALPX= RAO*ALP 
GD TC lAO 
10 WRITE (N2»180) 

Ntw = l 

ALP » 100. 

CL = 100. 

C ♦♦♦♦NON-ANSI**** 

READ (N5.P) 

LN * RN*l.E-t>*.t 
TXT « 3HALP 

IF (KOOE.bO.O) TXT - 3H CL 
CALL SECJND(TIME) 

WRITE (N2. 200) TXT. CL A i flODb +i ) » LN» N, NS» T IME. Kb L0» P.Cl» kuE L/ 

1 RSCP.SETA.ST.PCH.SEPM. XSEP.NPTS. IS.LL. 12 

2 .EPSl.NOES.REX. NFAST.NPELAX.ITYP.FSYH.TSTEP.DEP 
IF (ABS(XOUT) .CT..;>) 6J 10 7727 

GO TO 7726 

7727 CALL OUTPT 

CALL PL0T(0..0..999) 

STOP 

7728 CONTINUE 

C SELECT OUTPUT TAPE 

N<r B 

IF (lZ.Gb.80) NA s N2 
C2 » .5*(GAMMA-1 . ) 

C7 = GAMMA/(GAMb'A-l. ) 

IF ( ALP.EO.IOO. ) GO TO 20 
C ALP HAS BtEN INPUTTbL. KEEP IT FIXED 

NCY = 0 
MODE « 1 
ALPX » ALP 
20 ALP = ALPX/RaD 

IF (CL. EC. 100.) GO TC 25 
C CL HAS BEEN INPUTTED, KEEP IT FIXED 

NCY « 0 
MODE = 0 

YA = .5*CL/CH0-DPHI 
DO 114 L * 1,M 
DO 114 J - ),NN 

114 PHI(L,J) = PHKL, J)+YA*PhIP(L) 

DPHI » .5*CL/CHL 
CLX = CL 
26 CL » CLX 

C CHANGE PAkAHETEPS whlCH DE^^END ON IHb MACH NUMBER 

EM = AMAXKbM, .1C-4C) 

IF (EM.NE.EHX) NCY = 0 
Cl » C2+1./(EM *£M ) 

C6 = C2+EM *EM 
C4 = 1.+C6 
C5 = l./(C0*C7) 



QCftIT » (C1+CI)/(GA«HA+1. ) 

BET » S0RT(1.-Eh *EM )-l. 

C CHECK FOK TERMINATE# K£ TRiEVt^ DR STOKE INSTRUCTIONS 
C IK WILL BE -1 ONLY IF THERE IS A NAMILIST ERROR 

IF ((ITYP.EG.OJ.GR.( IK. EQ. -in GO TO 1?U 
CALL COSI 

IF (KS.NE.O) bO TO AC 
REWIND N3 

IF ( ITYP.GT.O) GO TO 30 

WRITE! N3) COHC# PHI# AA# bB» AR COLD# ANGOLO# XOl 0» YUL D# UE LULO# k# K S # R I 
1 #0SUM# GAMMA, XMON# R3C P# RFLC# R DE L# BC P# NS 1# KP# ST 
GO TO lAO 

30 READ (N3) COMC # PHI# AA# BB# ARCO LD# ANGCILD# XQLD# YOL D# DE LOLD# k» kS# R I 
1 #DSUM#GAMHA#XMCN#RBCP#RFL0#RDEL#BCP#NS1#KP#ST 
CALL MAP 
GO TO lAO 
AO CONTINUE 

IF (NS.GT.O) 60 TO 70 


70 


99 

100 


NS 

GU 

IF 

60 

IF 

GO 


80 


lOS 


lAl 


1A7 

IBl 


TO CRUDE GRID IF ITYP.GT.O 
(ITYP.GT.O) CALL kEMESH(-l) 

TO lAO 

(ITYP.GT.O GO TO 99 
BACK TO FINER GRID 
CALL REMESH(l) 

GO TU lAO 

SET UP CONSTANTS AND 03 CUNFORMAL MAPPING 
KD = 1 
XPHII >- 0. 

IF ( RCL.NE.O.) XPHll = 2.+CHU/RCL 
XA » 1.-2./RFL0 
AN60 = -RAO^BBd) 

TXT = 3H CL 

IF (MODE.tC.O) TXT = 3HALP 
DO AT MOS' NS CYCLFS 
IF (RN.Lh.O.) NSl * IOjOCOO 
IXX s M+2 
IXX » lXX-1 

IF (XC(IXX-I).GT.XMON) Gj TO 8C 
LC » 0 

DO 12C K • 1,NS 

IF (NDES.GE.O) 60 TO lOS 

IF (MOD(LC»S6).NE.O) GJ T3 lOS 

WRITE (N2#210) TXT 

LC « LC+1 

CONTINUE 

IF(NFAST.LE.O) 60 TO 
CALL SWEEPl 
IF(NRELAX.LE.O) GO TO 
DC 1A2 Lr«l#NRELAX 
CALL SWEEP 
CONTINUE 
NbW=*0 

NCY = NCY+1 


lAl 

1>1 




^SPuoii' 
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ALPX - RAOiALP 
CLX« 2.+OPHI+CHD 
YA « YA*XPHII 

C WRITE REidOUALS ON EVERY KP CYCLES 

If (N0ES.6E.0) GO TO UO 
IF (MOO(K»KF) ,NE.0) 60 TO 110 
LC « LC+l 
IN0C0*1 

CALL GTUPB(D1»D2,CP1»C0W#SL»RDEL»PUCP) 

IN0CC«0 

ANGO » -RA0*3B(1) 

WRITE (N2>190) NC Y# YR/ YA# 01, C2# IK# JK, NSP»CL A ( 2-MUOE ) / ANGC> CP1» 

1 CDW#CD4 

C DO A 60UNDRY LAYER CCRRECTION EVERY NSl CYCLES 

110 IF {nOD(K#NSl).NE.O> 60 TO 12i» 

IF (K.EQ.NS) GO TO 140 
WRITE (N2#190) 

LC « LC+1 
FSYH « 6. 

CALL 6TUKE(D1#02#CP1# tJCP#SL#ROEL#RbCP) 

ANGO » -RA0*BB(1) 

IF (hODE.EO.O) OPHI = .5*CLX/CHD 
C CHECK TO SEE IF wE HAVE SATISFIED CONVERGENCl CRITERIA 

125 IF ( AMAXl ( ABS (YF )#ABS(YA)) .LT.ST) GO TC 310 
120 CONTINUE 

310 IF(NDES.LE.O) 60 TO 140 
CALL CYCLE 

IF (M00( K0#K0ES I.RE.O) GO TO 138 
Call GTURB(01#02#CP1#6CP#SL#F OEL#KtCP ) 

CALL rtAP 
CALL COSI 

138 IF (KD.EO.NOFi) GC TO 139 
KD a KD+1 

GC TO 100 

139 REWIND N3 
ITYP = 1 

WRITE (N3) COMC#PHI# AA,BB# ARCOLD# ANGOLDfXJLO# YOLD#OELOLD#R#RS#RI 
1 #0SUK#6AK'U# XMON#RBCP#RFLG#ROEL#dCP#NSl#KP#ST 

140 HYP a lAbS(lTYP) 

CL » CLX 

LN a RN*l.E-6+,5 

XPF a XPF*AMIN0(1#IABS(M4-N4) ) 

XP a XP*XPF 
CALL SECONOniP.E) 

NTPE a n 4 
TXT a 3HALP 

IF (HODt.EO.O) TXT = 3H CL 

150 WRITE (NTPE#200) EN# TXT# CL A ( MODE + 1 ) # LN# rt#N# NS# T INE # FF LO# t« C L » M)E L# 

1 kBCP#PcTA#ST#PCH#S£P«#XSEP#NPTi# IS#LL# 12 

2 #EPSl#N0ES#kEH#NFA31#NR£LAX#lTYP#FSYM#riltP#ULP 
IF (NTPE.E0.N2) GU TO l-sO 

NTPE a N2 
GO TO 150 

160 IF (ITYP.GE.2) CALL GTOkd ( Dl# D2# CPI# BCP# SL . Ki)E L# R bCP ) 
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tMX • EM 

nyp«i 
6Q TO 10 
170 ITYP » 4 

IF (iK.EC.-i) WRITE (N^>220) 

CALL <>TURB(01 *02»CP1# dCPi iL / Kt)EL#Rl3CP ) 
terminate plot 

CmLL PLCT( 0«»0.»999) 

CALL EXIT 

180 FORMAT { 7H REAL' P/l 

190 PCRMATI t>X> I9» <»F12.i# 1««^ I J, I (3/ c FI 1 «3. Fil 51 

,3X5HR6^S:Ff? P P LCF . . 2 , 3 X. 9HRCL« F 4 . 2, 3X.HKL't L - f 9 . 3 

3xJHlZ-I3/^4H^tD5?^Ju''^‘^^’^'^^^^^^’^*^^*^*^P'‘5«12»3X3riLL.13, 
3X3H1Z-I3/ 6H bPSl*Fd.4»3X6H NLES«I3/3X 

,4HkEM»F8.4.7H NFAST»I3,3X,7HNRELAX = I3W, 

-.1 l>Hf S Y.i«F3 . 4»3Xf6HrSTfcP»>f6.4»3X»4HDtP«.FH u. / i\ 

JJO fD««»I <21H0...N.nELUr t«PQE...,10r,2CHE(UGR4i; lu lEKKINirt ) 


1 

2 

3 

4 
b 
b 


SUBRUUTlNk LSTERR 

CuMMON /A/ M(47)»IK 

IK » -1 

RETURN 

END 


bUBR'JUTlNE RESTRT 

COMMON PHI tl62»31).FP(I62»31)^A(31)»tj(3H#r(31)>0(31)/E(31) 

\ AA( 162), 86(162 )>C0(lt2) 

COMMON /A/ PI»TF,RAO,CM,ALP,KN,PCH,XP,K,CHD,OPHI,CL,RCL,yK 
1 ^XA, YA, TE»DT, OR, DlL Til/ DLL R,.v A, DON, 1)5 N> RA4,£ PS I L,LCR IT, Cl, C/; 

Jptf t, r . -CO I; ""■'‘'■‘^'‘'^‘^^^‘''^^''''"^121M,N2,N3,N4,NT,1xX 

,NPTS,LL,I,LjEP,H4,NLw,EPS1,NDES,aLEN,jCAL0i 

, sc A LOO, No, GAMMA, NOPT, CSTAk, REM,i)£P,01NF,TSTEP,XE,uT 
,INC,OFAC,GAM,KDE:>,PLIiZ,OPL,OPO *iiltl,XLUl 

SET UP CJNSTANTS 
TP = PI+PI 
KAO * 180. /PI 
ALP = ALP/RAO 

IF ( ( N ♦ I ) , N E • NN , UR . ( M + 1 ) . N t . M M ) NCY = C 
MM = M+1 
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IF (LL.E0*0) LL ■ M/2+1 
NN ■ N+1 

DK » -l./FLOAT(N) 

OT » TP/FLOAT (K) 

OCN « COi(OT) 

DSN « SlN(OT) 

DELS « .5/OR 

OELTh • ,5/DT 

RA ■ OT/DR 

RAA * DT*DT 

00 10 K « 1>N 

R(K) • l.+OR+FLDAT(K-l) 

RS(KI » (RA*R(K) )+(RA*R(K) I 
RHK) * -.25*0T/R(K) 

10 CONTINUE 
R(NN) « 0. 

BtT ■ SORT( l.-EM+EM) -1. 

IF (NOES.GE.O) 60 TO 5 
CALL AIRFQL 
GO TO 6 

5 CALL REAOQS 

6 CONTINUE 

IF (MOOE.EQ.l) CL « 6. +P 1*CHD*S I ( 1) / ( 1 , +3P T ) 
OPHI • .5+CL/CHO 
HA « MM/3 

Mb « HH-2+ ( (MA+1 ) /2 ) 
IF((NT.GT.1A0).0R.(XP.LT.0.)) JK - -1 
J«1 

00 AO L » 1#HM 
OELOLOIL) » 0. 

OSUH(J) » 0. 

ARCOLO(L')»ARCL(J) 

IF(J.GE.MM) 60 TO 70 
IF((J«LE«MA). OR .(J.Gb.iiB)) J*J + 1 
OSUH(J) « 0. 

J»J + 1 

AO CONTINUE 
70 NT ■ L 


WRITE (NA,100) NT 

100 FORMAT (1H0,IA,A5H POINTS WILL BE USED TO DEFlNt INNER AlRFOll 
InTPL(NT^ARCOLD» X0L0/AKCL/XC» PHI(l/i)/PMl(l,5)>PHl(i,7) ) 
INTPL(NT/APCULD» YOLUf ARCL#YC»Pril(l,3).PHI(i,5),PHI(i,7) ) 

SPLIF{HK,ARCL,FM,PH1(1,3),PHI(1,5),PH1(1W),Lc.,Lo.) 
lNTPL{NT»ARCOLD#ANGOLD>ARCLiFH,PHI(i, J),PHI(l,5),PHi(l,7) ) 

L ■ 1» M 


50 

60 


CALL 
CALL 
CALL 
CALL 
CALL 
00 60 


DO 50 J 

PHI(L.J) 

CONTINUE 


1/ NN 

‘ R( J)*CO( L)+i)PriI + PHlR( L) 


rSYM 
IS • 2 
RETURN 
END 


FSYM-12, 
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subroutine cosi 

SINES>C0SINES/ ANO THE TERM AT INFINITY 

01311, t(jii 

1 ^ ^ ^ ^ 31 ) > R ( 31 ) ^ RS ( 31 ) # R 1 ( 31 ) ^ AA ( 162 ) > BB ( 1ft? I rnfi^?i 

COMMON /A/ PI ^ R AO# EM# ALP# KN# PCH/ XP#TC#CHO# DP hT • r I on 

1 >*?'J*>Tfc><'T,DI<,OUIH,OtLI<,H,I,CN,osS,RM“Ejsu!ScRn^t; ?? 

• ItJlM, N2>N3>NA#NT# I XX 

>NPTS,LL,I,LSEP,HA,NEw,£PS1,NDES,XLEN,SCALQI 
>SCAL00,N6»GAMMA,NQPT,CSTAR, REM,UE?,Q1NF,TSTEP,X0UT 
.lNC.OFAC,GAM,KDfS,PLTSZ,OPL,GPU ^ EP» XDUT 

TPI = l./TP 
ANG « ALP + OBd) 

SN s SIN(ANG) 

CN = SORT (l.-SN+SN) 

DO 10 L « 1,H 
COIL) « CN 
SKL) = SN 

SN a CO(L)*OSN+SN*OCN 
ANG a ANG+OT 
10 CONTINUE 
CO (MM) a CN 
C0(MM+1) a C0(2) 

SI(MM) a SN 
SKMM + l) a si(2) 

RETURN 

END 


SUBROUTINE SWEEP 

SWEEP THROUGH THE GRIi) ONE TIME 

COMMON /FL /FLU XT-<» » CD4 


COMMON PHI(162»31),FP(l62,31),A(31),B{31),CI31).Dt^i » 

1 »RP«31»>RPP(31),R(31)>RS(31),PI(31),A1(J62 r? 

2 'Sl(i62),PHIR(l52),XC(io2),YC(162),Frt'({i2)’AR?uJt2dDSUMf 

COMMON /A/ PI» TP» RAD#EM> ALP#RN# PCri/ XP» TC» CHO# DPHI/CL»RCL.Yk 
^ *^*'][^'^^'®^'*^’’'*^^*-TP'^DELR/RA/CCN# 0 SN»RA<V/EPSIL/QCRIT Ci r? 

.C^,C5,C6,C7,BET,BEIA,FSYM,XSEP,SEPM,dLE(A),MU%.NN 

>IK»JK»IZ^ITYP^M0DE»lS»NFC^NCY»N><iN»NG»IDIH»N2iN'‘? Mt k.t*t 

>NPTS,LLd,LSEP,M.,NEW.EPSl,ND[s,uSiCArd^ 

>SCAL00#N6»GAMMA/N0pT>CSTAR^REtl^0EP»QINF»TSTFP.YniiT 
.INC,0FAC,GAM,KDES,PLTSZ,QPL,QPU 


YP 


0 . 


NSP a 0 

00 10 J a i,nN 



» PHI(1»J)+0PHI 
PHI(r.H<-l#J) a PH1(2#JJ+0PHI 
J ) « 0» 

RP4(J) a 0. 

RP5(J) « 0. 

10 RPP(J) a 0. 

SWEEP THROUGH THE GRID FROM NOSE TO TAIL ON UPPER SURFACE 
T c ■ “2 • 

LLP»U+1 
DO 30 I«LLP>H 
CALL MURHAN 
DO 30 J a l,N 

30 PHI(I-1,J) a PHU I-l, J)-Rp (J) 

DO 32 J>1^N 

32 PHI{M,J)«PhI(M,J)-E(J) 

DO 51 Jai,N 
E(J)aO. 

RPP( J)«0, 

RP4(J) a 0. 

RP5(J) a 0. 

51 CONTINUE 

SWEEP THROUGH THE Gkli) FROM NOSE TO TAIL UN LOWER SURFACE 

• t — C 0 

I » LL 
80 I a I-l 

CALL MURMAN 
DO 60 J a 

60 PHI(1+I,J) a ?Hin+l,J)-RP(J ) 

IF (I.GT.2) GO TO 80 
DO 70 J a 

70 PHI(2,J) = PHI(2> J)-£(J ) 

OC 11 J- IfNN 

PHI(MM+l»J)aPHI(2»J)+DPHI 

EU)aO. 

RPA(J) a 0. 

KP5(J) a 0. 

11 CONTINUE 
TEa-2. 

I«HH 

CALL MURHAN 
DO 50 Jal,N 

PHI(MM>J)aPHI(HM»J)-E(J) 

50 PHKl, J)apHl(MH» J)- 0 PH 1 
DO 12 Jal,N 
b(J)aO. 

RP4(J) a 0, 

PP5(J) a 0. 

12 CONTINUE 
TEa2. 

I=LL 

CALL MURHAN 
Du 13 Jai^N 

13 PHI{LL/J)aPHI(Ll,J)-E{j) 

ADJUST CIRCULATiON TC SATISFY THE KUTTA CUNOITIQN 
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Ml 

42 

90 

95 

97 


242 


100 


If <HooE!eia,'‘JS*}o‘5J‘‘2'ii‘OPHi)).ciEirH.st(iii 
Jf' *^I^^5.6E,0) go to 41 
ALP « ALP-.5+YA 
60 TO 42 

BB(1) . BB(l)-.5*yA 
CALL COSI 
GO TO 95 

YA ■ TP*YA/( I.+bET) 

DPHI * OPHI-t-YA 
DO 97 L « 1,M 

FLUXT^S!” ’ '’'’“^♦PHIRCL) 

NF«N-10 

IP(N*LT«30) NF«N*5 
DO 242 L*2«HN 

QF»(U*U+V*V)/FPCL,NF) <L»NF-1 n*0£LR -CO(U 

FLUX.RH*V/R(NF) 

fluxt-fluxt+flux 

CONTINUE 

fluxt»dt*fluxt*cho 

FLUXT4«FLUXT 

JJ*(*°de.eq.0) return 
00 100 J » 1,N 
00 100 L • 1,M 

Sr^'^’ ’ PHI(L,J)«YA*PHIR(U 
ENO 


c 

c 


c 


JUDKuuliNE nURMAN 

3 ^ANG0LD(162)>X0LD( 162)^ 

^ ^RP4(31),RP5(31) ‘^‘^^^’"^'^C°L0(162),0EL0L0(162) 

^Z,ITYP,MOOE,IS,NFCj.NCY»NRN 

if 

E(NN) « 0. 

FAC ■ -,5*TE 


113 



BU 

QS 

cs 


IH « I-l 

IM ■ I+i 

KK s 0 

PHIO * PHK I# 2 )- 2 ,* 0 K*C 0 ( I ) 

PHIYP« PHI(I, 2 )“PH 1 (I, 1 } 

PHIYY » PHIYP+PHIO-PHI ( I, 1 ) 

PHIXP « PHI(l+l, 2 )-PhI(I-l, 2 ) 

CHECK FOR the tail POINT 
IF (I.NE.MH) 6 C TO 10 
C(l) * (C 1 +C 1 )*RS( 1 ) 

A(l) = -C( 1 )+XA*C 1 -C 1 

10 U « PHlXM*OfcLTH“SI(I) 

' U/FP(I,1) 

• U*BC 
' C1-C2*QS 

BO - aQ*OS*(FP(]-l,l)-^p{I+l,l,, 

X » RA 4 *(CS+ 0 S)*CD(I) 

C(l) = (CS*CS)*PS( 1 ) 

* CS*kS( 1 )*PHIYY + KI( 1 )*hCi + x 

CMOS ■ CS-QS 

PHIXT = ttt TA* ABS ( U) +ABS ( CdOS ) 

EM 2 s QS/CS 

EPS 2 = EPSl*VLAYEP(th 2 » 0 PL,UPU) 

PHIXXX. EPS2*PH1XX-RP4(1) 

RPA(l) = EPS2*PHIXX 

0 ( 1 ) = 0 ( 1 ) + PhIXXX 

IF (wS.Lfc.OCPIT) GU TJ 30 

FLOW IS SUPERSONiC, BACKwAwD LlFFbKtNCES 

PHIXT * PhIXT-CHQS 
PHIXXH = pop(l) 

A(l) = “(C ( 1 ) ♦PHIXT) 

Ed) - 0(l)+CHCS*PHIxXM-PHlXT*r(l) 

A(l) = A(l)- 2 .*fcPS 2 -RP;;(i) 

“psin 

GO TO <»0 

*n subcritical/ central differences 

•jO A( 1 ) a XA^CMOS “C(I)-PHIXT 

D(l) a D( 1 ) +CM 0 S*PHIXX-PHIXT*E ( 1 ) 

Ad) » Ad )- 2 ,«EPS 2 -RP 5 (l) 

0(1) -( EPj2+2 . *RP3 ( 1) ) ( 1 ) 

= EPS2 

OCi non-boonoaky points 

^0 RPP(l) a PHIXX 
00 60 J a 2 >N 

PHIXY a PHIXP-PhiXM+(c(Jn)-L(J-l))*FAC 


Ad) 
0 ( 1 ) 
RPb(l) 
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45 


PHIXH « DU 
DU « DU*U6LTH 
PHIYYM a PHIYY 
PHIYH » PHIYP 

PHIYP « PHI(I,j*l)-Phl{i,j ) 
PhlYY B PHIYP-FHIYf« 

U « R(J)*&U-SI(I» 

RAV » R(J)*RA*V 
BO » l./FP(I,J) 

UCU « BO^^U 


LS 

UV 

VS 

cs 


BQU«U 

{ BOU-f BUU) *V 
BQ*V*V 
US + VS 
CS = C1-C2+US 
C^^VS « CS-VS 
CMuS a CS-US 
PHIXT a BETA*A3S(U) 

PHIYT = BE1A«ABS(PAV» 

EM2 « OS/CS 

EPS2 = EPSi♦VLAYER(t-^^2,QPL,QPb) 

PHIXXX= EPS2*PhIXX-PP^( J ) 

RP4(J) = EPS2+PHIXX 

1 FP I + l,J))+PAVM»-.p(l,j-i,.Fp(i,j^ ,,, 

D(J» = D(J» +PHIXXX 
UV a .5*BCU*RAV 

G3 TO 50 

USE BACKWAkD DlFFtRENCING 


IF (OS.LE.CCRin 
SUP6RSUN1C FLU// 
KK a KK+1 
CMOS a CS-OS 

FO a 1 ./ 3 S 

AUU 
BUU 

evv 

AVV 
BUV 
AUV 


US + FO 
PS ( J )*AUU 
VS*FO 

ps{j)*evv 

UV*^FO 

BCJ«ABS(RAV)’f'fC*rE 
» BVV*PHlXX-BUV*PHIXY+tUU*PHIYY 

‘ C 5 *BUU 

a PHIXr-ChOS*(AUU+AUu-AUV) +CS*BVV 
a PHIYT -CMOi ♦ ( A VV>A VV-AUV ) 

C(J) a 3 (J)+PHIYT 
PHIXXM a PPP(J) 

IF (V.LF.O) GO TO 43 

PhpYM a PHI(I,J+ 2 )-PHI(l,j+i)_phXYP 

PHIXYM a PHIYP+Hhl ( Ih, ( IM, J+l ) 

GU TO 46 

BfJ) a C(J) 


PHINN 
B(J) > 
PHIX I 
PhlYT 


TLP .15 

1 /J)- 
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C(J) > 60 

46 PHISS ■ AtU^PHlXXM-fAUV^PHI aYM-»AVV*PH 1 YYM 
A(J) » -(b( J)*C(J)+PhlXT) 

D(J) « 0( JJ ♦CMiS*PHlSS+CS*PHINN-t(J)*PHlXT 
A(J) « A(J)-2.*EPS2-KP5{J) 

D(J) • 0(J) “(tPS2+2.*<P5( J) )*E( J) 

RP5(J) » EPS2 
6D TU 60 

C SUBSONIC FLOw» LSE CENTRAL uIFEEPENCtS 
50 C{JI = RS(J)*CFVS 
B(J) ■ C(J)*PHIYT 
PHIXT « PHIXT+CMUS 
A(J) « XA*CKUS-B(J)-C(J)-PH1XT 

D(J) » 0(J)+CMUS*PHlXX“UV*PhlXY+C {J )*PHIYY“PHIXT*E( J) 
A(J) ■ A( J)-2.*EPS2-KP5( J) 

DU) » DU) -(EFS2+2.*KP5(J))t=l (J) 

RF5(J) « EPS2 

IF ( V.LT.O.) GO TO 60 

fi(J) « CU) 

C(J) ■ C(J)+PH1YT 
60 RPP(J) » PHIXX 
NSP « NSP+KK 

C SOLVE THE TRIC1A60NAL jYSTEN 

CALL TRIO 
RETURN 
END 


C 


C 


C 


SUBROUTINE TRIO 


SOLVE N DIMENSIONAL TRIOIAGONaL SYSTEM OF EQUATIONS 
COMMCN Phl(lb2/31)»FP(162/ 31)>A(31)»B(31)/C(Jl)/i}(31)/E(il) 

1 /RP(31)>KPP(3I)#M3l)»RS(31)>Rl(3l)/AA(lc2)>6B(l62)#C3{lo2J 

2 >SI(162)#PHIR(le2)»XC(lb2)»YC(162)>FM(162)»AftCL(162)»CSUM(lo2) 

3 » ANG0LD(162) #XCL0(162)/ YJLD(152)>ARCGLD(lb2),0EL0LU{l62 ) 

4 >RP^(3I)»RP5(31) 

COMMON /A/ PI/Tr/RAD#EN>ALP#RNiPCri>XP»TC»CHD»OPHI>CLiRCL/YR 

1 >XA>YA»TE»DT>0P»UELTH/DELR#PA»DCN»0SN>kA4»EPS1L>UCR1T>C1»C2 

2 »C4,C5»C6,C7»4ET»BETA»FSYM#XSLP,SePM»TTLE(4)>M,n,NM,NN,,nSP 

3 > IK> JK» I2/ITYP»M0UE» IS>NFC>NCY# NKN» NG> IJIM»N2/ N3»N4>NT# IXX 

4 >NPTS»LL#I/LSEP# M4» NEV*/EPSl>N0LS>XLbN>SCAL01 

5 >SCAL0a>N6#GAMhA>N0PT#CSTAP> FFM/DEP>01NF# TSTEP» XOUT 

6 t INC»0FAC»3AM,KDES»PLTSZ»aPL>QPU 

XX = 1./A(1) 

PP(1) « Ed) 

Ed) » XX*D(1) 

DO ELIMINATION 


DO 10 J » 2>N 

C(J-l) » C(J-1)*XX 

XX * l./(AU)-B(J)*C(J-l) ) 

RP(J) » E(J) 

10 E(J) = (D(J)-B(J)*b(j-i))*XX 
DO BACK SUBSTITUTION 


Reproducibility of the 
SS PA«E IS POOR 
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EMX « ABS(E(N}) 

00 20 J * 2*N 
L « NN-J 

E(LJ « E(U-C(U*E(L+1) 

20 EMX » AMAXKEMX# ABS(E(L) ) > 

FIND THE LOCATION OF THE MAXIMUM RESIDUAL 
IF (EMX.LE.ABS(YR) ) RETURN 
IK > I 

DO 70 J » 1,N 

IF ( ABS(E( J)) .EQ«EHX) 60 TO 74 
70 CONTINUE 
74 JK » J 
YR - E(JK) 

RETURN 

END 


C 

C 


SUBROUTINE REMESH(LSIGN) 

GO TO CRUDER GRID IF LSIGN IS -1 
GO TO FINER GRID IF LSIGN IS *1 

COMMON PHI(162»31)»FP(162»31)>A(31)> B(’li.C(31)/D(3l)>E(31) 

9 AA(162 )» BB ( 162 )> CO (lb2 ) 

COMMON /A/ PI/TPaRAD>EM»ALP>RN#PCH#XP»TC»CHOaDPHI,CL>RCL> YR 
^ 'Z^'^^'°^'‘^‘*'^E‘-TH>DELRaRA,DCN>0SN»RA4,EPSIL>QCRITjC1>C2 

ISaNFC#NCY,NRNaNGaIDIH»N2,N3#N4»NT#IXX 
>NPTS# LL> I» LSEP>M4»N£Wy EPS1#NDES>XLEN»SCALQI 
^SCALOOa N6»GAHHAaNQPT» CST AR»REHaDEP>QINF» TSTEP»X0UT 
/ INC#QFACaGAMaKDESa PLTSZaQPLaQPU 


X » 2. ♦♦LSIGN 

N6 » FL0AT{N6)/X+.2 

M « FLOAT(M)^X +.2 

N « FL0AT(N)^X+.2 

IF (N*E0*14) N>15 

LL « FL0AT{LL-1)^X+1.2 

IF (LSIGN. GT.O) MM » M+1 

IF (LSIGN. GT.O) NN » N+1 

LSEP » FL0AT(LSEP-1)^X+1,2 

PF » 1,/X 

DELR « X^DELR 

DELTH = X^DELTH 

DR • PF*OR 

DT « PF^OT 


) 


OCN « COS(DT) 
DSN » SIN(DT) 
RA4 • PF^PF^RA4 
NCY » 0 
I » LSIGN 
MP » HM+1 


F 
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CALL PERMUT (R»NN>1) 

CALL PERHUT (RS,NN,1) 

DO 5 J B 1>N 
5 RI(J) ■ -.2i*DT/R(J) 

CALL PtRHUT (OSLM,MP/l) 

00 2C L B If UN 

20 CALL PERMUT ( PHI (1# L )/MP, 1 ) 

DO 30 L * 1>MP 

30 CALL PERMUT ( PHI { L> 1 )# NN/ IDIM) 

MM B 

NN B N+l 

IF (X.EC..5) GO TO HO 
00 40 L B 1»M»2 

DSUH(L*1) * .5*(CSUM(L»+0SUM( L+2) ) 

DO 40 J « 1/NN,2 

40 PhI(L + l»J) B .5*(PHI(L»J JtPHI(L + 2» J) ) 
DO 50 J * 1,N,2 
DO 50 L » 1»HM 

50 PhI(L»J+l) = .5*(PHI(L»J)+PHI(L#J+2)> 
80 CALL MAP 
RETURN 
END 


SUBROUTINE PEPMUT (AX,NX#JX) 

REORDERS POINTS wITHiN AN ARRAY 

COMMON PH I (162* 31 )> FP ( lo2» 31 ) > A ( 31 ) » t (31 )» C ( 31 ) » D ( 31 ) , b ( il ) 

1 /RP(31)/kPP(31i/P(31),RS{31),kI(31), AA(r62),Bb(i62)/CC(102) 

2 >SI(162)»PHlR{l62)#XC(162)>YC(lb2)/FMU62)#ARCL(lo2)»USuh(lft2) 

3 # ANGOLO( 162)# XCLD(lo2)» YOLO (162) » ARC OLD (162 )/ DELClD(1o2) 

4 »RP4(31)#RP5(31) 

COMMON /A/ PI#TP#RaO#E. 1#ALP#PN# PCH# XP#TC#UhU»UPHl#CL/RCL# Yk 

1 #XA#YA>Tb#DT#OR#CbLTH#OELK#fiA#DCN# 3SN#KA4#tPSIL#wCRIT#Cl/C2 

2 »C4»C5#C6#C7#3t.T .PbTA#FSYM» XS c P# S fc PM# T T Lb ( 4 ) # M# N, MH> NN# N SP 

3 # IK# JK# 12# ITYP#NOLE» 1S»NFC»NCY» NkN» NG# I DI M , N2» N3# .N4» NT# i X X 

4 #NPTS#LL# I#LSEP#M4,NEb#EPS 1#NLES#XLEN# iCALOl 

5 # SCAL0Ci#N 6#GAMMA#NQPT»CSTAP# PEM# Dt P# uiNF# T S T LP# XLUT 

6 # INC#QFAC#GAM#KDbS#PLrS2»UPL#0PU 
DIMENSION AX(1) 

L B 1 

JY B Jx+jx 

NY B 2*((NX-l)/2)+l 

NZ B 2*(NX/2) 

IHI.GT.O) GO TO 30 
NY B JX*(NY-1)+1 
NZ B JX*(NZ-1) 

DO 10 J = 1»NY#JY 
A(L) B AX(J) 

10 L « L*1 

DO 2C J B JX#NZ#JY 
A(L) * AX(J+1) 

20 L B L+1 
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6U TO 60 

30 OCi J = 1»NY»2 
A(J) s AX(U 
40 L » L^JX 

OU 50 J « 2»NZ>2 
A(J) < AX(L) 

50 L » L+JX 
60 L » 1 

00 7C J » 1#NX 
AX(L) a A(J» 

70 L « L+JX 
RETURN 
END 


SUBROUTINE GETCPICOF) 

C COMPUTE CP»CD» AND CM oY INTEGkATIuN AND OUTPUT MACH OIAGkAN 

COMION/FL/hLUXTA^CDA^COw# INUCD 

COMMON PHI(16?#31)»FP(lo2»31)#A(3i)»3(31)#C(Jl)#D(31)»K3i) 

1 * RP (31 RPP( 31 )> k( 31 )» RS( 31 ) >RI { 3i )> AA( io2 )» BB( 162 )#C0(io2 > 

2 /SI(io2)»PHIP(162)»XC(162)iYC(ib2)»FM(l62)» APCL(l62)»DSoN(lfc2) 

3 f ANGOLO( 162J » XuLL (162 )»yJL0(l62)>AKC0LD(162)^D£LCLD(lo2) 

4 ,P?A(31)#RP5(31) 

COMMON /A/ PI»TF/RAD,E 1,ALF#RN>PCH, XP#TC/ChJ,DPHl/CL>RCL# Y k 

1 »XAiYA>TE»DT/DR#D£LTH»0£LR/PA#DCNfl)iN/<A^/EPSlL>QCKlI/Ci»C2 

2 #C4>C5»C6»C7>BE IxJETAf FiYM#XSEF>SLPM;TTLE(n)>'1,.,,M.'1,Ni\^NSP 

3 > IK>JK» IZ» ITYP»MU0E> IS»NFC»NCY»N*»h^ N6» IDIM»N2»N3>N4>NT» IaX 

4 /NPTS/LL#I/LS£P>H4f N£„ tPSl,NLbS/ XLcN^SCALQl 

5 } SC A LOO# No# GAMMA# NOP r# CS TAR# RfcM# UbP# wlNF# Ifcp, xuUT 

6 #INC#aFAC#GAM#KLES#PLlS2#0PL#UPU 
RLAL MACHN»MACh 

COMPLEX CLCD#TMP 

DIMENSION MACHn(1)#CPX(1)#MN(1) # 1M^CH(21) 

EOUI VALENCE ( MACHN ( 1 ) # A ( 1 ) )#(CPX(l)#PHIRli))#(MN(l)#FP(l#ji)) 

DATA IM ACri/lHC# IHft# lH3#i.Hr#lHL# IHV# IriR# IHX# IHY# IHZ# IhO# Inl# lh2# 1H3 
1# 1H4#1H5#1H6# lH7#lH0#lH9#lri+/ 

Data tx /4hcdf=/ 

MACH(G) » SaRT(C/(Cl-C2+U ) 

IMC(C) = MIN0(21# IPIX(1U.»C)+1) 

CLCO a 0. 

CM a 0. 

IP ( ( XP.GT.O. ) .CR.( I2.LE .30) ) GO Tu 10 
OY a YOLO(NT)-Yr>LD(l) 

IF (FSYM.NE.O.) DY = YC (MM )-YC ( 1 ) 

REWIND M4 

wRITb (M4#12C) LM#CL#UY»TC> NRN»MM 
10 LO 20 L a 1#MM 
CP a CPX(L ) 

C COMPUTE CP+02 

TMP a CP<=SJkT(F?( L» 1 ) ) +CMFLX( CU3( PM( L) ) » S1N( FM (L) ) ) 

C SUM UP CL# CO# AND CM 

CLCO a CLCO+TMP 


WEPRODUCIBILITY op the 
ORIGINAL PAGE IS POOR 
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c 


c 


c 


CM = CH-KXC (L )-.25 )*ftEAL ( IMP }-YC { L ) *AIKA3( TMP) 
WHTE PUNCH OUTPUT UN .i'r Ih XP«0 AND IZ.mT« 60 
IP ({XP.GT.0.).0fi.(I2.LE.80)) GQ TO 20 
Q ■ MACHN<L)*SOkT(Cl/(l,+C2»^ AChN(L)*KACKN(L) J) 
V • 0*$1N(FM(U) 

U » 0*C0S(FM(L)» 

IF (XP.EU.O) GO TO lt> 

WkITE (M4/130) U»V»XC(Ll^yClD(L)>CP 
60 TO 20 


15 bRITE (MA#130) U# V> XC { L)> YC ( L CP 
20 CONTINUE 


CORRECT CL»CD FOR ANGLE OF ATTACK 

CLCD a -(0 T*CHD)*CLCU*CMPLX(S1N(ALP)>C05(ALP)) 

CM » OT*ChD*CM 

WRITE CD»CL/CM ONTO N4 

COW * REAL(CLCD) 

QCR»SORT(QCRIT) 

0CD4=2.*{UCR-1.)*FLUXTG 
CD « CDW+COf 
CD4aC0+0CD4 
CDWaCOW+OCOR 


IF ( INOCD. ECitO ) PRINT 291# COU# CD F> CuR 
261 F0RMAT(5H CDta»F10 . t># 5H CDF «F10. 5# <.H C0»F10.5) 
CL2 » AIMAG(CLCO> 

IF (INDCD.EQ.O) GO TO 16u 

CALL CDSI 

RETURN 

160 CONTINUE 

IF (MR.EW.N3) GO TO 

IF { CDF. EO.O. ) GO TO 70 

WRITE (N4#90) EM# CL2# CM# COw» T X# CDF# CDR 

GO TO 80 

70 WRITE (N4#90) £M#CL2#CM#CDA 
C CONSTRUCT MACH NUMBER DIAGRAM 

WRITE (NA#1A0) 

80 I » IMC(EM) 

I - IMACH(I) 

C USE PRINT WIDTH OF 12 FOR MACH NUMaER DIAGRAM 

KB a MM 


MC a MAX0( 1#M6/IZ) 

HA a MC+MaXO( 1#MB-IZ*MC) 

C WRITE OUT MACH NUMBERS AT INFINITY 

WRITE (N4#100) (I# L a HA#Mb#MC) 

C DO MACH NUMBERS ONE LINE AT A TIKE OOrtN TO ThE BODY 

J a NN-MC 

AO RSJ a R(J)tp(J) 

DO 50 L a MA#Mb#KC 

U a (PHI(L+1#J)-PHI(L-1#J) )*R( J)f'OELTH-SI(L) 

V a (PHI(L#J+l>-PhI(L# J-1) )«OELP*RSJ -CU(u) 

0 a (UfU+V^V) /FF(L# J) 

1 a IMC(MACH(Q)) 

MN(L) a IMACH(I) 

50 CONTINUE 

WRITE (NA#100) (HN(L)#L * MA,MB#HC) 
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J • J-MC 

IF (J.6T.1) GO TO <»0 
C DO TH£ LINE WHICH :S THE 40DY 

00 6C L > HA>l16#nC 

1 « 1MC(HACHM(L)» 

60 MN(L) » INACH(i) 

WRITE (NA»100) (NN(L)>L = «A,Hh>h':) 

IF (ITYP.GE.A) CALL GRAFIC(CD) 

RETURN 

85 RNX « .l*AINT(RN*l.E-5) 

WRITE (NA#150) £H>CL> TC»CH» RNX>Cl)F 
RETURN 

90 FORMAT ( 1H12X iHEM = F5 .A » AX3HCL«F7 . AX3hCM=F6 . A» AXAHCOW=F 7 . A A AA 
1 »F7.5»AX iHC0»F7.5/// ) 

100 FORMAT (3X/13CA1) 

120 FORMAT (3H M»> FA , 3» 5X» JHCL=» F5 # 3> 5X» 3HC Y=# F5 . 3/ 6X» AHT/C=> 

1 FA.3»lAXf 215) 

130 FORMAT (A020) 
lAO FOR'IAT (IhO//) 

150 FORMAT { IHO// 7X3hEMa» FA ,3» AX3HCL«> F6 . A# AXAHT /C»# FA ,3> AX3HCM*# 

1 F6.A#AX3HRN*»EA.l, AXAHCUF»»F6.A/) 

END 


SUBROUTINE GRAFIC(CD) 

COMPLEX ZP>ZO#SFAC»SIG 
REAL MACHN 

COMMON/FL/FLUXTA,CDA#CDw» IN CCD 

COMMON PHI(162#31)»FP(162>31)^A(31)>fa(31)>C(31)»i){31)>E(31) 

1 >RP(31)»RPP(31)>R(31)»PS(31)>RI(31)>AA(i62)/BB(162)*C0(ic2) 

2 »SI (162)»PHIR(lb2)#XC(162)#YC(lfa?)>FM(i62)»ARCL(it2)/DSUM(le2) 

3 » ANG0LD(162)»X0L()(l62)>Y0LD(lO2)> ARCOLO( 162)^DEL0L0{ 162) 

A >RPA(31 )»RP5 (31 ) 

Common /a/ PI»TP>RA[,»E1, ALP/RN» PCH, XP,TC/CHO»OPHl>CL#RCLf Yk 

1 f XA/ YA» Th»OT#OR#OELTH»OELR/RA^ DCN# DSN» R AA# E PS I L » OC R I T»C 1 » C2 

2 f CA>L5#f 6>C7»6Er»bETA> FSYM,XSEP>SEPM#TTLE (A)>M,N#MM,NN,nSP 

3 » IK^ JK> IZ» ITYP»MCDE> li»NFC»NCY»NRN» N6» 1 01 N2> N 3> N A» N ( » IXX 
A >NPTSf LL>I/L3EP#MA»Nb*.» EPSl#Nr.bS>XLbN»SCALQl 

5 t SCAL0''*N6#GAMMA#NUPT»CSTAK» KLM/DbP/ ol 'lF>TSTEP#xCUT 

6 #INC» C/GArt,KOES>PLrS2»OPL>OPU 
DIMENSION CPX(l),MACHN(l)>T(o) 

EQUIVALENCE ( CPX ( 1 ) , PH IR ( 1 ) ) #( MACHN ( 1 )» A( 1 ) ) 

DATA TOL/l.b-o/ $ PF/-.A/ » SC F / 5 . 0/ » YCR / A. 0// S I Zb / . 1 A/> SCO/ toO . / 
C MOVE THE ORIGIN TwO INCHES OVER AND TWC INChcS UP 

CALL PL0T(2.0>2.5»“3) 

Y0R = AMAXi(3.5t .5*AINT(20. ♦EM-7.0) ) 

C PLOT CP CURVE AS A FUNCTION Lib X 

CPF » l./PF 
CCP » CPF*CPX(1) 

C«LL PL0T(SCF*XC(l)»Y0P+CCP/3) 

DO 10 L = 2»MM 

Cup « AMIN1(6.5-YCR»CPF«C0X{L) ) 
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10 CALL PL0T(SCF*XC(L)» Y0K+CCP#2) 

C DRAW ANi) LABEL ThE CP-AXiS 

CALL CPAXIS{-.5#YCP» l.-l./PF»7*b~Y0R»PF ) 

C CUMPUTE AND PLCT CRITICAL SPEEC 

CALL SYhbOL ( -. YOR ♦CPF‘»C PX ( MM + 1 ) / 2 . *S IZE » 0. »-l ) 

C PLOT 0OOY 

CALL PLQT(SCF*XC(l)#SCF*YC(l>/3) 

DO 20 L - 2>PM 

20 CALL PL0T(SCF«XC(L)>5CF»YC(L)>2) 

C LABEL THE PLOT 

ALPX « RAO*ALP 
TXT=>6HANALY3IS 

IF (EPSl.GT.O.) TXT » lOHART. VISC. 

IF(FSYM.Gb.6. ) TXT=6HTHE0RY 
XL»-.9 

C ♦♦♦♦NON-ANSI - SEE WRITEUP AT END**^* 

IF(FSYM.GE.6. ) GO TO 30 

IF ( (NOES. LT.O) .AND. (EPSl.LE.O. ) ) GO TC 200 

TTLE(l) ■» AHVISC 

TTLE(2) « AHQUS 

TTLEO) * AHCESI 

TTLE(A) « AHGN 

ENCODE (oC#21C»T) TTLL»N»N>NCY> EPSl 
GO TO AO 

200 ENCOPE (60»191»T) TTLL»M#N»NCY 
GO TO AO 

30 LN«RN*l.E-o*. 5 

ENCODE (60> 190»T) TTLE»M#N»NCY/LN 
AO CALL SYK3t:L(-l.lA»-1.0#SIZE>T#O.»b6) 

C ♦♦♦♦NON-ANSI - SEE WRITEUP AT END**** 

ENCODE ( 60/ 170> T) T XT> EM> AL PXt CL f CDA 
IMCDA.LT.O) SNCCDE(fcOi 171» T) T XT» EM# AL PX/ CL» CDA 
CALL SYMBDL(XL#-1.3i#SI2E#T#C.#60) 

CALL SYMBQL{XL-.10>-1.35+.t?*SIZE»l.t)*SlZt#l5#0.»-l) 
CN=C0(1) 

SN=S1 (1 ) 

C READ AND PLOT EXPERIMENTAL DATA IF XP IS NOT ZERO 

IF (XP.Ew.O.) GO TO 130 
REWIND MA 
READ (MA#1A0) NP 
IF (EOF(MA).NE.O) GO TO 130 
IF (NDES.GE.O) GO TO 220 
READ (?;A»150) EhX# ALPX#CLX#CDX# SNX 
READ (MA,160) ( CO ( L ) # S I ( L » # L = 1»NP) 

TXT » lOHEXPERIMENT 
GO TO 230 

220 READ (MA#2A0) TCX# DCAVE# YRX, SNX 

READ (MA»160) ( CD( L ) # S I U ) # L « 1»NPJ 

TXT « 6HINPUT CP 
230 CONTINUE 
NC = b9 

1F(SNX.GE.0.)6G TO 50 

TXT»tHDESI6N OhRTWI; ?V\G£ IS 

NC«3 ^rw.Ujii7 
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c 


c 

c 


c 


c 


c 


c 

c 


♦♦♦♦NON-ANSI - SEE WHIfcUP AT END**^ + 

50 ENCODE (60»170*T) TXT# tHX» AL PX> CLX# CDX 

Tc **‘^*®’ CNCODF(bO# T1#T) T X T# EM X# ALPX# C LX# COX 
IF(NDkS»6E*0) ENCODE ( b0#<»50# T ) TXT# TC X# OOAVE# YRX 
CALL SY,1dCL(XL#-l .7»S12E»T#0.#60) 

CALL SYMBOL ( XL-. 10# -1, 7 ♦.i ♦SIZE# SIZE# NC#0.#-1J 

DO leO L • 1#NP 
CCP » YOR*CPF*SI(L) 

IF (CCP.6T.8.4) GO TO 180 

CALL SYMBOL (SCF*CO(L)» CCP# .5*S1ZE»NC#0.#-1) 

160 CONTINUE 

130 IF (1TYP.E0.5) 60 TO 122 
PLOT THE SONIC LINE 
EX « l.-EPSIL 

SET SINES AND COSINES FOR USE IN FOURIER SERIES 
MX • M/2 


com - 1, 

Sim - 0 . 

Du 60 L » 1#MX 

C0(L+1) » CO(L)*DCN-SI(L)*OSN 
CO(MH-U « C0(L+1) 

SKL + IJ » CO<L)*OSN+S1(L)*OCN 
60 SKiM-U » -SKL + l) 

DO 120 L * 2»M 

LOOK FOR SONIC POINTS ON THE BODY 
IF (MACHN(L).LT.l. » GO TO IIC 
IF (MACHNd-ll.GE.l.) GO TO 80 
IPEN « 3 

COMPUTE Z AT SONIC LINE ON BODY 
70 Kl s (MAChN(L ) -1, >/ ( MAuHNC L )-HACHN( L-1) ) 

ZP » CMPLX(XC(L)4R1^(XC(L-1»-XC{L))»YC{L)+R1+(YC(L-1) 
CALL PL0T(SCF*PEAL(ZP)#SCF*AIMA6(ZP)#IPEN) 

IF (IPEN.E0.2) GO TO 12u 
FIND THE SONIC LINE ALONG A RAY 
80 0 = MACHN(L) 

SX » SI(LJ*CN+SN*CO(L) 

CX a CO(L)*CN-SN^SI(L) 

FAC » .5+DR 

ZO « CMPLX(XC(L)#YC(L») 

DO 90 J » 1#N 
ZP a SFAC 
kj » R(J) 

QS a Q 

IF (J.EO.l) 60 TO 82 

U a (PHKL + l# J)-PHI(L-1#J) )fRJ«OELTH-SX 
V a ( PHI ( L# J*1 )-PHI ( L# J-1 ) )«DELk*RJ*RJ-CX 
0 a (U^U + V*V) /FP (L» J ) 

0 a SQRT(U/(C1-C2*U) ) 

82 SIG a CMPLX(RJ^CO(L )#RJ^SI (L) ) 

COMPUTE ((1-SI6MA)++(1-EPSIL) )SIGMA 

SFAC a CEXP(EX*CL0G((1,#0.)-S1G))/SI6 

SUM UP FOURIER SERIES TO OBTAIN CONJUGATE OF N 

S a -Bed) 

00 8A K a 1,NFC 


YC ( L ) ) ) 
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c 

c 

c 


84 

86 


90 

100 


LT ■ H00((L-1)*K»M) 

Sj'_S*RJ*{AA(K+l)*SI(LT+l)-BB(K*l)»CO(LTnn 

IF (RJ.LT.TOL) GO TO 86 
CONTINUE 

COMPUTE THE ARGUMENT OF DZ/OR 

SFAC • -SFAC*CMPLX(COS(S)»SIN{S))/CABS(SFAC) 

magnitude to obtain OZ/OR 
SFAC « SFAC*(CHD<'SQRT(FP(L> J) ) )/(R( J)*R( J) ) 

PERFORM THE INTEGRATION 
ZQ ■ ZQ+FAC*SFAC 
FAC - OR 

IF (Q.LE.l.) GO TO 100 
CONTINUE 


110 

120 

122 


ZQ > 
ZP • 
R1 • 
ZP « 
CALL 
GO TO 


Z0-,5*0R*SFAC 

ZQ-.5*DR*(SFAC+ZP) 

(Q-l.)/(Q-OS) 

ZQ*R1«{ZP-Z0) 

PLOT (SCF*REAL(2P»#AMAX1(-2.0>SCF*AIMAG(ZP) )»2) 
120 


SURFACE DELS #0.xl9) 


132 


134 


140 

150 

160 

170 

171 

190 

191 
240 
210 


IPEN o 2 

IF (HACHN(L-l»*Gt.l«| GO TO 70 
CONTINUE 

POSITION PEN AT BEGINNING OF NEXT PAGE 
CALL FRAME 

CALL PL0T(-2.0,-2.5>-3) 

IF <{FSYH.NE.7.).0R.(ITYP.EQ.6)) RETUft., 

PLOT THE BOUNDARY LAYER DISPLACEMENT 
INOEXR (0.»XC>H) 

PL0T(2,»1.5»-3) 

SYMB0L(1.36>-.65>SIZE/19HL0WER 
CPAXIS (0.#0./0.»4.^ l./SCD) 

LOWER SURFACE 

PLOT (SCF*XC(l)>iCD*DSUM(l)#3) 

2»MX 

(SCF*XC(L)>SCD-t-0SUM(L)#2) 

PL0T(0.»4.5.«-3) 

SYMBOUl.36,-,65^ SIZE/19HUPPER 
CPAXIS (0.>0.#0.#4.»1,/SCD) 

UPPER SURFACE 

PLOT (SCF*XC{MX)>SCD*DSUM(MX),3) 

DO 134 L - MXfH 

CALL PLOT (SCF*XC(L+l)/SC0*0SUM(L+l)/2» 

CALL PL0T(10.>-6.>-3) 

RETURN 

(10X,I3) 

(3F6.3#F7.5,E9.1» 

(2F10.4) 

( A12#4H M»F4,3#3X4HALP«F5.2» 3X3HCL«# F*i.3.3X3HCD«. Fb.4) 

^•*'''^*3,3X4HALP-F5.2,3X3HCL-,F5.3,2X3HCD-F6.4)^^ 

MILLION) 

cnnu*^*'**'*'^’‘^^”*^“^^'^*^*^2'3X4HNCY«14»4X12HN0 VISCOSITY) 
FORMAT (F7.3#2E10.2>F4.1) 

FDRHAT(4A4,3X4HM+N*13,1H*I2,3X4HNCY«I4,4X,5HEPS1»/F5.3) 


MX 

CALL 
CALL 
CALL 
PLOT 
CALL 
DO 132 L ■ 
CALL PLOT 
CALL 
CALL 
CALL 
PLOT 
CALL 


FORMAT 

FORMAT 

FORMAT 

FORMAT 


SURFACE DELS >0.»19) 


124 



ri o o o r> o 


250 F0RMAT(A12»2X#AHT/C«F5.3»2X#3HD0»E8.2»2X»5H0PHI«E8.2» 
END 


SUBROUTINE CPAXIS ( XOR# YOR» BOT/ TOP> SCF ) 

DRAWS ANO LABELS THE CP A/IS 

XOR»YOR IS THE LOCATION OF THE ORIGIN OF THE AXIS 
BOT IS THE LENGTH OF THE AXIS BELOW THE ORIGIN 
SCF IS A SCALE FACTOR USED FOR LABELING 
DRAW THE LINE FOR THE AXIS 

SCF NEGATIVE FOR CP AXIS ANO POSITIVE FOR DELS AXIS 
SIZE • .I2-SI6N(.02»SCF) 

CALL PLOT (X0R»Y0R+T0P#3) 

CALL PLOT (X0R»Y0R-B0T»2) 

C DRAW HATCH HARKS ANO LABELS ONE INCH APART 

N ■ 1+INT(B0T)*INT(T0P) 

S * -AINT(BOT)*SCF ♦I.E-IZ 
XH « X0R-(3,*SI2E)/.7 
YH « YOR-AINT (BOT) 

00 10 I > 1*N 

CALL SYMBOL ( XOR# YH, SI2E» 15y 0. #-l ) 

C ♦♦**NON~ANSI - SEE WRITEUP AT THE END**** 

IF (SCF.GT.O*) ENCODE (10,25, A) S 
IF (SCF.LE.O.) ENCODE (10,20,A) S 
S ■ S+SCF 

CALL SYMBOL ( XH, YH, SIZE, A, 0., A ) 

10 YH « YH+1. 

IF (SCF.GT.O.) GO TO 30 

CALL SYHB0L(X0R+.l,Y0R+2.5, .1A,1HC,0.,1) 

CALL SYHBOL(XOR+.25,YOR+2.38, .14,1HP,0.,1) 

RETURN 

C DRAW THE X-AXIS 

30 CALL PLOT ( XOR, YOR-BOT, 3 ) 

CALL PLOT (XOR+5.0, YOR-BOT, 2) 

CALL SYMBOL ( XOR + 5. 5, Y0R-.07, , 1<», IHX, 0., 1 ) 

YH ■ Y0R-B0T-SIZE-SI2E 
DO 40 I ■ 1,5 
S « .2*FL0AT(I) 

ENCODE (10, 20, A) S 

XH » Y0R<-FL0AT(I)-SI2E-S1ZE 

CALL SYMBOL ( XH, YH, SIZE, A, 0., 4) 

40 CALL SYMBOL ( XOR+FLOAT ( I ), YOR-BOT, S I ZE, 15,90., -1 ) 
CALL SYMBOL ( XQR+.25, YOR+3 .0, . 14, 4HDE LS, 0. , 4 ) 

RETURN 

25 FORMAT ( F4.3) 

20 FORMAT (F4.1) 

END 
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SUBROUTINE 60PL0T (NRN) 

INITIATE PLOT 

***********^:^**^***0^^^^,*^^t*******^** *4 ********* ************* 

THIS SUBROUTINE SHOULD BE REPLACED BY ANY ROUTINE WHICH INSTRUCTS 
THE SYSTEM TO INITIATE A PLOT 

***************************************************************** 
IF (NRN. 6T. 1000) 60 TO 50 

CALL PL0TS(600#26HJEFF MCFADDEN# FOLDER 3210) 

RETURN 

50 CALL PL0TSBL(600#26HJEFF MCFADDEN# FOLCER 3210) 

RETURN 

END 


C 


1 

2 

3 

A 

1 

2 

3 

A 

5 

6 

1 

1 

2 

3 

A 


SUBROUTINE AIRFOL 

READS IN DATA FOR AIRFOIL AND MAKES INITIAL GUESS FDR MAPPiNO 
FUNCTION BY COMPUTING FOURIER COEFFICIENTS 
IF ONLY X>Y COORDINATES ARE PRESCRIBED SMOOTHING IS DONE 
COMMON PHI(162>31)#FP(162/31)>A{ 31/# b(31)»C( Jl)>0(3i)»E(3i) 
#RP(31)#RP?( 31)#R(31)»RS{31 )»RI(31) » AA ( 162 ) # B6 ( 162 ) # CO ( lfa2 ) 

»SI (162),PHIR(162)»XC(162)#YC(162)»FM(162)#ARCL(162)#DSU«(162) 
# ANG0L0(162) »XOLD(lb2)#YaLD(l62)#ARCOLO(162)#DELDLD( 162) 
#RPA<31)#RP5(31) 

COMMON /A/ PI#TP#PAD#tM#ALP#kN#PCH#XP#TC»CHO#DPHI,CL#RCL#YK 
»XA»YA#TE#0T»0r.,DELTH,0ELR#RA»DCN#DSN»RAA,EPSlL»tCRIT»Cl#C2 


#CA#C5#C6#C7#8ET#DETA#FSYM#XSEP#SEPM#TTLE(AJ,H#N#MM#NN#NiP 
#IK# JK»IZ#ITYP#HDDE# IS# NFC# NCY # NRN# NG# 10IM#N?#N3 NA#NT#1XX 
#NPTS#LL# I#LSEP#MA#NEW#EPSl#NDES#XLbN#SCALOI 
#SCAL30»N6#6AMMA»N&PT»CSrAR> REM# DEP# OINF# TSTLP# XLUT 
# INC#QFAC#GAH#KDES#PLTSZ#UPL#OPU 


DIMENSION XX(1)#YY(1)#U(D# V(D# 
#DS(1)#SS(1)#CX(1)#SX(1)#QSP(1) 
EQUIVALENCfc(XX(l)»FP(l#3 ) )#(YY( 
<V(l)#FP(l#7))#(#i{l)#rP(l#9) )# ( 
(1# 13) )» (TH( 1),FP(1#15) )# (TT(1 ) 
(SS(1)#FP(1#21) )# (CX(i)#FP(l#23 
FP(1#27))#(Z(1)#FP(1,29)) 

S0(Q2) = C2*Q2 

SM0QTH(Q1»02#03#QA) « 02+SC(SL(S 
DlSCOl) = (Q1-ERR)»( (U1-ERR)’«=(01 
DATA T0L#NT#ISYM#C0NST»VAL/.AE-7 
DATA OXOS1#OXDS2#OYDS1#DYOS2/A+0 
NmP is the NUMBEk OF POINTS IN C 
LC » NFC 
NMP « 2*LC 
MC = NMP + 1 
PILC « P1/FL0AT(LC) 

IF (FSYM.GE.6.) GO TO 150 
WRITE (NA#A70) 

REWIND N3 

READ «N3»A10) TITLE 
IF (FSYM.GE.3.) GO TO 100 


W(i 

)# SP( l)#CiRCi 

(1) 

#TH(D# 

TT(1) 

#TI 

TLE(15)#Z(1) 





1)# 

FP(1,5) )# (UC 

L )# 

FP( 

i#l) 

}s 

SP( 

1)#FP(1#11) )i 

.(L 

IRC 

(D# 

hp 

# FP 

(1#17))#(DSC 

D# 

FP( 

1#19 

))f 

) }# 

(SX(1)#FP(1#25) 

)#( 

QSR( 

1)9 

O('JA) ) )*.25*(01- 

-02 

-C2 

■••U3) 


-ERR)+CONST) 





#999# 0# .2#AriRUN 

/ 




./ 

# XT/-1./ 





IRC 

LE PLANt FDR 

FOURI 

ER S/RIES 
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^0 

c 

50 

100 

110 

120 

130 

140 

95 

150 

160 

170 


re«o^in coo«oinmes as PKauuiEf BY PSUGSAns 0 AND f 

XXd) » o! 

NL » 2 
rewind N3 

READ (N3>510) EM, Ct> DY, TC, NKN 
IMC s HOOi INT (100,#CM + .5 >, 100 ) 

ICLl a M00( lNT(CL+.05), 10) 

ICL2 » M00( INT (10 .*CL+. 5 ) , 10) 

ITCl a MJOnNT(10.^TC + .05),10) 

1TC2 a MO0(INT(10O.*TC+.5),10) 

MOOE^^ «^0*i»30,mt) lflC,ICLl,lCL2,ITCl,ITC2 

IF (NRN.LT.O) FSYM«2. 

DO 40 L a 1,999 

read (N3,500) U( L ), V ( L ) , XX ( L ) , YY( L ) » F AC 
♦♦♦♦CHECK FOP END uF FlLE«**t t‘«-J,i-AC 
IF ( E0F(N3) .Ht.O) GO TO 50 
IF ( XX(L ) .LT. XX(NL) ) NL a I 

continue ) 

AlRFfjIL HAS BEEN EXTENDED IN PROGRAM O 

NT a L_i 

NRN a IaBS(NFN) 

60 TO 150 

READ IN AIRFOIL DATA FkOil CARDS 
read (N3,420) FNU,FNL,EPSIL 
read (N3,470) 

NT a FNU+FNL-1. 

NL a fNL 

DO 110 1 a nL,NT 

read (N3,4?0) 0(1 ),V{I),XX(I),YY(I) 
read (N3,470) 

DO 120 I a 1 ,nl 
J a NL+l-I 

U(J),V(J),XX(J),YY(J) 

Lu 130 J « Ifif 
TTLE(J) a riTLE(J) 

IF (FSYM.LE,4.) GO TO 150 
DO 140 L a 1,NT 
TH(L) a XX(L)/kAO 
XX(L) a u(L) 

YY(L) a V(L) 

GO TO 195 

epsil*^^^o stfeam function 

If'^hLym^fo^i^^ lengths Can EE computed to FIRS 

UU 160 I a 1,NT 
Th(I) » 0. 

ISYM a I 
Gu TO 200 

COMPUTE SLOPES FROM VELOCITIES 
TH( 1) a ATAN( V(1 ) /U(l ) ) 

OSR(l) a 0 ( 1) #U( 1 ) +V ( 1 ) ♦V ( 1 ) 


CiROER 
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OQ 190 I • 2>NT 

C CHOOSE NEAREST BRANCH EOR THE ARCTANGENT 

TH?I) + I)) > 

190 QSR(i) » u(n*un)+vm*v(i) 

195 IP (EPSILtCT.l.) EPSIL ■ ( TH( 1) -( P I*TH( NT ) ) ) /PI 

IF (FSYM.GT.5.) EPSU « ( TH ( 1 ) ♦ TH( 2) -Th{NT )-TH(NT-l ) ) /TP-1. 

C COMPUTE ARC LENGTH TO FOURTH ORDER ACCURACY 

200 spri) « 0 . 

DO 210 I - 2>NT 

OUM ■ AHAXl(.lE-20».5*A8S(TH(I)-TH(I-l)n 
OX » XXdJ-XXd-l) 

DY « YY{I)-YY(I-1) 

210 SP(I) ■ S?(I-l)+SORT(DX*OX+OY*OY)*uUH/SiN(OUM) 

ARC » SP(NT) 

SN « 2. /ARC 
SCALE « .25*ARC 
EE ■ «5*(1. -EPSIL) 

DO 220 L » 1»NT 
220 SS(L) » ACOSd.-SN*SP(L) ) 

SS(NT) a PI 

IF (ISYH.NE.O) GO TO 350 

CALL SPLIF (NT>SSdH,U,V>W,3»0.,3>0.) 

IF {FSYM.6T.5.) GO TO 232 
WRITE (NA#410) TITLE, VAL>NRN 
IF (N4.NE.N2) WRITE (N2,A10) TI TLE, VAL, NRN 
C PRINT OUT DATA ON THE AIRFOIL 

WRITE (NA,A30) 

DO 230 L a l/NT 
VAL « TH(L)*RAD 

SUM »-SN*U(L)/AMAXl( .1E-5,SIN{SS(L) ) ) 

IF ( (L *E0.1) .OR* (L.EC.NT) ) SUM « V( L ) *S 1 GN ( SN, FLOAT ( L-2 ) ) 

230 WRITE (NA,A80) XX ( L), Y Y (L ) , SP (U , VAL, SLM, V ( L ), W (L ) 

WRITE (NA,AA0) 

C MAKE INITIAL GUESS OF ARC LENGTH AS A FUNCTION OF CIRCLE ANGLE 

232 DX a (XX(NT)-XX( 1) )/TP 
DY a (YY(NT)-YY(l) )/TP 
DO 2A0 I a i,MC 
AN6L a FLOATt X-1)*PILC 
CIRC(I) a ANGL 
CX(I) a COS(ANGL) 

SX(H a SIN(ANGL) 

YY(I) a 1. 

IF (EE.NE.O.) YY(I) « (2.-2.*CX ( I ) ) ♦♦EE 
FAC a SIGN(1.+CX(I),FL0AT(LC-I) ) 

2A0 SP(I) a AC0S(.5^FAC) 

3P{HC) a PI 
CIRC(MC) a TP 
IF (FSYM.LT.O.) GO fU 2AA 
SCALE a ARC/ARCL(MM) 

SNlaZ./ ARCL(MM) 

DO 322 I«1,H 

322 ARCL ( I )aACOS ( 1.-SN1+ARCL( 1 ) ) 

ARCL (MM) api 
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250 


260 


DUM 

CA 

DB 


00 2A2 L • 1,MM 
2A2 2(L) « FL0AT( L-1 )*0T 

2-^4 00 2<,r[ i 

BB(L) ■ CX(2*L-1J 
2^5 AA(L) • -SX<2*L-1) 

o“S tloT’ ITJoY'"""' ” M.FriCIENrs 

ensure Closure 

DUM ® 0« 

SUM ■ 0« 

FAC « 0. 

DO 260 L * 1»NHP 
DIM B OUH -TT(L) 

SUM « SUM-TT{L)*CX(L) 

FAC B FAC*TT(L)*SX(L) 

' 0U^1/FL0AT(N^^P) 

00 290 I « 1, NrtP 
SUfI = OS(I) 

290 DS(I) « YY{I»+EXP(SUM) 

OS(MC) = DS(l) 

Z ( 1 ) *0 . 

VAL* .5*PILC 
VALl»PILC/3. 

Z(2)=\/AL*(0S{1)+DS{2)) 

NiBNK + 1 
DO 295 J=3#NI#2 

Z(HC)«0. 

Z(MC-1)=VAL*{ DS(MC) +LS(MC-i ) ) 

NI I»NFC-2 
00 299 J*2#N1I,? 

HCJbNC-J 

Ho cSr 

ULi 3C1 J = 3^NI»2 

DS1 = Z(NFC+J)-Z(NFC*J-1.) 


295 

296 
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Z(NFC+J-l)»2tNFC4-J-2)-Zl 
IFtJ.EC.NI) 60 TO 303 
21«Z (NFC+J+D-ZtNFC+J) 

303 CONTINUE 

2{NFC+J)-Z(NFC+J-1)-DS1 
301 CONTINUE 

SCALE » ARC/Z(NC) 

ERR » 0. 

00 310 I « 1#NMP 

VAL = Acos(i.-2.=>2( n/z(MC n 
ERR = AMAXl(LRR#A3S(SP(I»-VALn 
310 SP(I) = VAL 

If (fSYH.LE.5.) WRITE (NA#A90) ERR/0A,0e 
IF (ERR.LT.TCL) GO TL 33o 
320 CONTINUE 

WRITE (N9/450) 

330 CALL FOUCF(NNP,TT»CX,iiiJ»AA) 

AA(1) - AkC 

AA(2) • l.-EPSIL-(0X*SlN(B3(l))+Dr*CGS(33(im/SCALE 
BB(2) » (-Dx*cos(etj(in+Dr*siNus(in )/scale 

IF (FSYM.GT.5.) GO TO 3A2 
WRITE (N4#460) EPSU» «MP 

IF ( (FSYM.NE.1.).ANC.(FSYN.NE.3.)) GO TO 3^1 
DO 3AA L » 1»KM 
3AA ZIL) » FL0AT(L-1)*0T 

CALL SPLIF(NC»ClKC»SP/U»Vi w^3>C.#3,0.) 

CALL TNTPL(HN,Z#DS/CIRC>bP#U> V»W) 

CALL SPLIF (NT»SS»QSk»U#V#w#1^0.#l»0. ) 

CALL INTPL(Kh»DS/»A#SS> 3SR>U» V# W) 

DO A L = i»MM 

A IF (A(L).LE.O.) A(L) - 0. 

3A1 IF (IZ.NE.120) GO TO 3A2 
•.RITE (NA,6A0) 

00 3A0 L » 1»NFC 
3A0 .iRITt (NA,A90) AA(L)>SJ(L» 

3A2 CALL MAP 
RETURN 

350 IF (FSYM.LE.5.) GO TC 355 
DXOSl » (XX{2)-Xx(in/3S<2) 

DXDS2 » (XX(NT)-XX{NT-U)/(SS(NT)-SS(NT-i)) 

OYOS' » (YY(2)-YY(in/5S(2) 

DYDS2 = (YY(NT)-YY(NT-1))/(SS(NT)-S3(NT-1M 
355 CALL SPLIF(MT»SS»XX#U,SP>W>l/DXOSlf l»DX0S2) 

CALL SPL1F(NT, SS» YY>V»TTiOS#l#DYDSif 1/DYDS2) 

IF (IS.LT.C) 60 TO 397 
OC * PI/FLUAT(NKP) 

ERR » SS(NL) 

DUN « DISCO.) 

FAC « PI/(DIS(PI)-UUH) 

DO 360 L » 1»MC 

360 CIRCCU » FAC*(01S(FLOAT{L-1)*DC)-DUM) 

CALL INTPL(NNP»CIPC»SX,SS^XX,U>SP^W) 

CALL lNTPL(r4HP»ClRC»CX»Si» YY/ V/TT#0S5 
SX(MC) = XX(NT) 
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c 


CX(MC) « YY(Nn 

SFAC » l^/iXXCND-XXCNL) ) 

XXNL » XX(NU 

DO 370 L « 1»MC 

CX(L) - SFAC*CX(U 

SX{L) » SFAC* (SM D-XXNL) 

XX(L) * SX(L) 

370 YY(L) * CX(L) 


wHITf (N^#p20) is 


IF (N2*NE«NA) WRITE (N2#530) IS 
IF(IS.EQ.O) Gu TU 39t> 

00 IS SfOOTHING ITERATIONS 
DC 3,90 K * 1, IS 
DO 380 L » 2»NMP 

XX(U » SMUOTH(SX(L-l),SX(L),SXiLtl),SX(U) 
380 YY(L) » SM00TH(CX{L-1),CX(U,CX(L+1),SXIU) 
00 390 L » 2»N«P 


SX(L) » XX(L; 
390 CX(L) s YY(U 
395 NT » HC 


CALL SPLIF(NT»ClPC*XX,U»SP»W»l,0.fl,G,} 
707 Tct!: ^'‘“‘'‘'^^*CIPC,YY,V,rT,DS,l,0.,I,0.) 

jVf liYM « 0 

IF (FSYM.GT.S.) GU TO 170 
U(l) = SP(1) 

V(l) = TT(I) 

U(NT) « SP(NT) 

VINT) = TT(NT) 

GC TO 170 


410 F0R.1AT (1X16AAW4) 
420 FORMAT (5F10.7) 


430 

44 0 
450 
460 


470 

4?0 

490 

500 

510 

520 

530 

540 


FbRMAT (35H0AIRF01L COORDINATES AND C OR ATUR t S /lhC» F-X. IhX . 1 n h y 
1 /9X#10HAFC LEN0TH,7X3hANG,bX5hKAPPA,lCX,2hKP,llX,3HRPP//) 

FORMAT (iHl»AX»3HEKl<‘»l«tX>2HDA>l^X>2riLa//) 

FORMAT (32H FOURIER SERIES DID NCT CONVERGE) 

INSIDE oF A CIkCLE//3X11H0Z/0SIoMa » 

1 jOFI 1/SIGMA**2 ) ♦( 1-5 IGMA) *♦( 1-EPSI L )*( LXP( W( S IGMA) ) // 3a» 

24cHw(SI3MA) s SL'h ( ' A(N)-1*B (N ) ) ♦SIGMA**IN-1 ) ) //3X# 7HsPSl L = 

3 F5 . 3» 20X> I4> 25F* POINTS AROUND THE CIRCLE ) 

FORMAT (IHl) 

format (F12.6»2F14.6/Fl^,3»F14t4#2E14.3) 
format (3E15.6) 

****CHANGE (4020) TO (20A4) UN IBM 3c0**** 

FORMAT (4020) 

FORMAT (3X»F4.3»dX»F5.3#8X»F5.3^10X,F4.i»l‘»X»I5) 

FORMAT (lOHOTHEi^E Art»IA,26h SMOOTHING ITEkAIIONS USED /) 
FGRMAT(4hAIPF»6X>3H01L/7X>I2>l>i->Il,bX/Il, 1H-# 211 ) 

FORMAT ( // 7X4HA(N)# ICXAHBI N) // ) 

END 


ce.- 


* 
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c 


c 

c 

c 


c 


c 


c 


SUBROUTINE NAP 

SUN UP FOURIER SERIES TO OBTAIN NAPPING FUNCTION 
CONPLEX TT#TNP _ 

CONNON PHI(162#31)#FP(162»31)#A(31)»Bf31)#C(31»#0(31J j-E(311l 

1 #RP{31»#RPP{31}/R(31)»RS(31)iRI(31)»AA(l62>#BBC162)#CO(H52) 

2 /SI(162)»PHIR(162)»XC(162)#YC(162»» FN(162)»ARCL(162)#0SUM(162 ) 

3 » ANGOLDI 162» >XOLO( 162 )» YOLO (162 )» ARC OLD (162 J»OELOLD( 162 ) 

A ,RPA(31)»RP5(31) 

CONNON /A/ PI,TP#RAD#EN»ALP»RN#PCH>XP>TC»CHO/OPHI#CL»RCL>YR 

1 ,XA,YA,T£#DT»OR>OELTH,OELR>RA#DCN,DSNj>RAA»EPSIL>OCRIT#Cl*C2 

2 ,CA#C5»C6/C7#BET>BETA»FSYNjXSEP»SEPM>TTLE( A)/N»N»HN»NN/NSP 

3 »IK»JK» IZ»ITYP»H00E»IS»NFC»NCY>NRN»NG»I0IN»N2>N3»N4/NT» IXX 
A #NPTS»LL#I#LSEP#NA/NEW»EPS1»NDES#XLEN»SCALQI 

5 ,SCAL00»N6»GANNA#NQPTf CSTAR#REN» DEP» QINF» TSTEP» XOUT 

6 ,INC#QFAC>6AN#KDES#PLTSZ>0PL>QPU 

♦♦♦♦CHANGE TO l.E-6 FOR SINGLE PRECISION IBN 360^^** 

DATA P0WfT0L/-12.»10.E-12/ ^ 

NOTE THAT THE SQUARE OF THE NAPPING NODULUS IS BEING CONFUTED 


NX • M/2 

SET THE SINES AND COSINES 
CO(l) > 1* 

SKI) • 0. 


DO 5 1 * 1»NX 

C0(L*1) “ CO(L)ADCN-SI(L)*OSN 
CO(NH-L) - C0(L+1) 

SKL + l) » C0( L)*OSN+SI(L)*OCN 
S SKHM-L) • -SKL-H) 

SET MAPPING MODULUS FOR CUSP AT THE TAIL 
DO 10 J • 1»N 

FP (1»J) * l.*R( J)*(R( J)-2.) 

DO 10 I = 1#MX 

10 FP(L*1»J) ■ l.+R(J)*(R( J)-2.*C0(L+1) ) 

IF (EPSIL.EO.O.) GO TO 30 

ADJUST IF THERE IS AN ANGLE AT THE TAIL 

DO 20 J * 1»N 

FP(1»J) “ FP(1#J)**(1»-EPSIL) 

DO 20 L - 1#MX 

20 FP(L+1»J) ■ FP(L*1>J)**(1»“EPSIL) 

NOW COMPUTE CONTRIBUTION FROM FOURIER SERIES 
30 DO 50 J « 1>N 

NFCX ■ HINO(NFC»1+INT(POW/ALDG10(R( J)-TOLI n 

RJ »• 2t*R(J) 

K « NFCX 


S » AA(K+1) 

35 S “ R(J)*S+AA(K) 

K » K-1 

IF (K.GT.l) GO TO 35 

FP(1»J) - FP(1» J)*EXP(S+RJ) 

DO 50 L - 1#MX 

K «• NFCX 

LX « K^L 

LT » M00(LX»M1 

S ■ AA(K+ )*C0(LT+1) 

Q « BB(K+1)+SI(LT+1) 
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40 LX « LX-L 

LT « M0D{LX>M) 

S * P(J)*S+AA(K)*C0(LT«-1) 

Q « R(J)*Ci + 8B(K)*SI(LT + l) 

K • K-1 

IF (K.6T.1) GO TO 40 
DUM » FP(L+1»J) 

FP{HK-L*J) « EXP(PJA(S-U) )*DUM 
50 FP(L+1>J) ■ FXP(PJ*(S+Q) )*OUH 
DO 65 L « 1#H 
S = PI-BB(l) 

DO 60 K » 1»NFC 
LT « M0D((L-1)*K>*1) 

60 S » S+AA(K+1)*SI(LT+1)-BB(K+1)*C0(LT+1) 

ANG = FL0aT(L-1)*DT 
FP(L#NN) » 1. 

65 FM(U « S“.5*(AK6+EPSIL*(ANG-PI ) ) 

FH{MP) » FH(l)-(l.+tPSIL)*PI 
DO 7C J = 1»NN 
FP{MH>J) * FP(i»J» 

70 FPt«M+l»J) » FP(2/J) 

C COMPUTE ARC LENGTH ANO BOOT FROM THE MAPPING BY INTEGPATluN 

XMIN « 0. 

YMIN * 0. 

YMAX = 0. 

S - -SQRr(FP(l»l) ) 

TUP * CMPLX(S*CCS<FM(1) )»S«SiN(FM(l) ) ) 

DO 80 L » 1»MM 
Q » SQRT(FP(L»1) ) 

S » S + Q 
ARCL(L) = S 
S = S+0 

TT = CMPLX(Q*COS(FM(L) )»Q*SIN(FM (L) » ) 

TMP » TMP+TT 
XC(L) » REAL(THP» 

YC(L) » A1MAG(TM?) 

XMIN » AMINKXMlN^REALaMP) J 
YMIN » AHINK YMIN/ AIMA3(TMP ) ) 

YMAX » AMAXKYMAX/AIKAGITMP)) 

TMP = TMP+TT 
80 CONTINUE 

CHD = l./( .5*XC{MH)-X'1IN) 

IC » ( YMAX-YMIN)*CH0 
DO ‘50 L = l/MH 
ARCL(L) = CHD*APCL(L) 

XC(L) = CHO*(XC(U-XMIN) 

90 YC(L) « CHD*YC{L) 

CHD = CHO/C.S^DT) 

IF tNDES.GE.C) RETURN 

IF ( ABS ( FSYM) .GT.5 . ) GJ TO IOC 

/NGO= -RAD*6B(1) 

WRITE (N4/120) TC/ANGO 

IF (N2.NE.N4) WRITE (N2/120) TC/ANGO 

IF (MODE.EQ.O) ALP = ( 1 , +Tt T ) *C L / ( 8 . * P I*ChO ) -BB ( 1 ) 
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100 CALL COSI 
RETURN 

120 FORMAT (32H0THE THICKNESS TO CHORD RATIO IS 
1 17H OF ZERO LIFT IS »F6.3»8H LEGRhES) 

END 


#F6.^//10H THE ANGL 


10 

20 

30 

40 


SUBROUTINE SPLIF ( N» S# F/ FP/ FPP# FPPP# KM# VM> KN, VN) 

‘ SUBROUTINE CONTRIBUTED BY ANTHONY JAMESON 
GIVEN S AND F AT N CORRESPONDING PCiINTS#COMPUTE A CUulL SPlINE 
THROUGH THESE POINTS SATISFYING AN END CONDITION IMPOSED On 
EITHER END. FP>FPP>FP?P WILL BE THE FIRST/SECONO AND THIRD 
DERIVATIVE RESPtClIVELY AT EACH POINT ON THE SPLINE 
KM IS THE DERIVATIVE IMPOSED AT THE START uF THE SFLINE 
VM WILL BE THE VALUE OF THE DERIVATIVE THEkE 
KN IS THE DERIVATIVE IMPOSED AT THE END OF THE SPLINE 
VN WILL BE THE VALUE OF THE DERIVATIVE THERE 
KMfKN CAN TAKE VALUES 1,2> OP 3 


MUNOTONIC 
S(l), F(l)/ 


S MUST BE 
DIMENSION 
K = 1 
M « 1 
I » M 
J » M+K 

DS » S(J)-S(I) 

0 « OS 

IF (DS.EO.O.) CALL ABORT 
DF » (F( J )~F( I) )/0S 
IF (IABS(KH)-2) 10»20,30 
U a .5 

V » 3.*(0F-VM)/DS 
GO TO 50 

U a 0. 

V a VM 
GO TO 50 


FP(i)» FPP(1)> FFPP(l) 


U 


- 1 . 


V a -OS*VM 
60 TO 50 

I a J 
J a J+K 

DS a SU)-S(I) 

IF (D+DS.LE.O.) CALL ABORT 
OF a {F(J)-F( I))/DS 
B a l./(OS+OS*U) 

U a 6*DS 

V a B<=(6.*DF-V) 

FPm a u 
FPpm * V 

U a (2.-U)*DS 

V a 6.*0F+0S*V 

IF (J.NE.N) GO TO 40 
IF (KN-2) 60#70»80 


fyit 
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frO V a (6.*v .-V)/U 
GO TO 90 
70 V a VN 
GU TO 90 

<30 V a f DS*VN+FPP( I ) ) / ( 1 , + pp ( 1 ) ) 

90 B a V 

P a OS 

100 OS a S(J)-S(I) 

U a f PP(1)-FP( I)+V 
FPPP(I) a (V-U)/DS 
PPPU ) a U 

V^a^U * ^ ^ ( V+L + 0) /6. 

J = I 
1 a I-K 

IF (J.NE.M) GO TO ICO 
FPPP(N) = FPPP(N-l) 

FPP(N) = B 

FP(N) . 0F+D*(FPP(N-1)+3 + 3j/(j. 

IF iKM.oT.O) RtTURN 

clmputl the imegkal in fppp 

V a FPP(J) 

105 1 a J 
J a J + K 

OS a S ( J )-S 'I ) 

U a FPP(J) 

FPPP(J) a FPPP(I)t.5*JSMF(l)+F(J)-0S*DSMU^V)/12.) 

IF (J.NE.N) GO TO 105 

RETURN 

End 


fil’PRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR 


10 


‘^*^»^I’f^*S,F,FP,hPP,FPPP, 

iFiiiiHsIs!" 

aS'lEMK!!:, 

J a C 

00 3C I a 1,NX 
VAL a C. 

SS a SKI) 

J a J + 1 


n a s(j)-ss 

IF (TT) i0/30>2C 
20 J a PiAXOUfJ-l) 
SS a SS-S(J) 


30 ?un 
RETURN 
HJO 
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SUBROUTINE CONJ ( N, F>6» X>CN, SN) 

CONJUGATION BY FAST FOURIER TRANSFORM 

GIVEN THE REAL PART F OF AN ANALYTIC FUNCTION UN THE UNIT CIRCLE 
THE IMAGINARY PART G IS C3NSTRUCTEU 
COMPLEX F»G»EIV,EIT 

DIMENSIUN F!1)>G(1)>X(1), CN(1),SN(1I 

DATA PI/3.lAlt><)26535b979/ 

L » N/2 

CX » l./FLOAT(L) 

ElV - CHPLX(COS(PI*DX)»SIN(PI*DX) ) 

DO 2 I » 1,L 
2 G(I) - F(1J 

CALL FFORH(L»G» X#CN» SN) 

G(l) * 0. 

1*1 

DO 1C J « 1»L>2 

EIT » CMPLX(SN( IX-OXjCNdjADX) 

1 « !♦! 

6(J) • G(J)*EIT 
10 G(J+1) « G{ J+1)*EIT*EIV 
DO 22 I«1>L 
22 SNd) « -SN(I) 

C*LL FTORM(L#G/X/CN/SN) 

DO 32 I«1»L 
32 SNd) » -SNd) 

EIV = CMPLX( AIMAG(G(L) )» REAL(Gd) > ) 

1 " L 

AO GCI) » Cf1PLX( AIMAG(Gd-l) )»REAL(Gd) ) ) 

I * I-l 

IF ( l.GT.l) GO TO AO 

G(l) s EIV 

RETURN 

END 


SUBROUTINE FOUCF ( N. G> X» A# B ) 

C FOURIER COEFFICIENTS BY FAST FOURIER TRANSFORM 

COMPLEX G,EIV,CP>X,GK 

DIMENSION G(l),Xd)> Ad)>Bd) 

DATA PI/3.1A159205358979/ 

L » N/2 

V « PI/L 

EIV » CMPLX(CCS(V),SIN(V) ) 

ENI ■ l./FLOAT(N) 

CALL FF0RM{L>6>X,A>B) 

GK » 0. 

I « 1 

DO 5 J = 1,L#2 

X(J) => CMPLX(B(I)»A(D) 

X(J + 1) = X(J)<»EIV 
5 I « 1 •♦•1 
K » L 
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00 10 J « 1,L 

OP » 6K-C0NJGI6(J)) 

GK » GK*CCNJG(6(J))-CP*X(J) 
ACJ) « -»(EAL(GK)*ENI 
b(J) • AI«AG(GK)*ENI 
GK « G(K) 

10 K * K-1 

A(L+1) » -B(l) 

B(l) . 0. 

B(L+1) . 0. 

RETUFN 

END 


C 

C 

C 


C 


5 

11 

10 

21 


SUBROUTINE FFORN( N, F, X, CN, SN ) 

FAST FOURIER TRANSFORM 

INPUT ARRAY F WITH REAL ANO IMAGINARY PARTS 
REPLACED BY ITS FOURIER TRANSFORM 
COMPLEX F(l),X{l),w 
DIMENSION CN(1)/SN(1) 

IF (N.LT.2) RETURN 
NS • 1 
NR « 2 
NO » N 


IN ALTERNATE CELLS 


SET THE SINES AND COSINES 
PI-3.141S92653S6979 
DT » (PI + Pn/fLOAKN) 

IF( ( SN( 1 ) ,£C.O. ) . AND. ( SN( 2 ) ,E0, S 1N( OT) > ) GO TO 11 


DO 5 J ■ 1»N 
CN(J) = COS(ANG) 
SN(J)*-S1N(ANG) 

ANG = ANG+DT 

DO 10 K « iMR, M 

If- (M0D(NQ#K) .ft.O) 60 TO 21 

CONTINUE 

ND*NO/K 


NS * NS*K 


NR ■ K 
IQ » 0 
ID » 0 

DO 22 I = 1,NS 
DO 2A J « 1»NC 
L » IQ+J 


LP=L+ND 

M=ID 

W «F(L) +F(LP) *CMPLX(CN(N+1),SN(H+1) ) 
lF(NR«£Qa2) 60 TO 2A 
L = LP 


DO 2b K»3»NP 
L ■ L+NO 
M a M + IO 
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IF (M*GF*N) M ■ M~N 
26 Vr » W+F(L)»CMPLX(CN(K*i),SN(M + l) ) 

2<t xno+j) = W 

ID » lO+NO 
1Q-I04-NC 

IFtlO.GE.N) IG»1Q-N 
22 CONTINUE 
NO > NO 

IF (ND.GF.l) 60 TO 61 
DO 32 K « 1>N 
32 F(K) » X(K) 

RETURN 

61 CO 60 K « NR»N 

IF (MOD(NO,K) .EC.O) GO TO 71 
60 CONTINUE 
71 NO=NQ/K 
NS » NS*K 
NR a K 
10 a 0 
ID a 0 

DO 72 I a i,NS 

DO 7<t J a 1,ND 
L a 10+J 
LPaL+ND 
M«10 

WaX( U+X(LP)*CMPLX{CN(M + 1), SN(N + 1 ) ) 

Ih(NP,E0.2) 6C TO 74 

L»LP 

DU 76 Kai,NR 
L a L+NO 

M a H+IO 

IF (H.GE.N) M a 

76 W a W + X(L)*CMPLX(CN{N + 1),SN{N + 1) ) 

74 F(IO+J) a w 
ID a ID+ND 
ICalO+NQ 

IF(IC.GE.N) IQaio-N 
72 CONTINUE 
NO a ND 

IF (ND.GT.l) 60 TO 11 

RETURN 

END 


Function inoexrix^apkay^n) 
dimension ARRAY(l) 

S a ABS(X-ARRAY(N) ) 

DO 10 L a l,^| 

IF ( ABS (X-ARRAY (L ) ) »UT. GO TO 10 
INDEXK a i 
S a ABS( X~ARRAY(L ) ) 

10 CONTINUE 
RETURN 
END 
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3 

4 
*> 
6 


10 


**^^'^8fDELPAX»DELBP/CPG>BCP/SL»R0EL#RBCPJ 
COMMON/ FL/FLUXT4»C04» COW/ INOCC 

COMMON PHI(l62/31).fP(l62/31)/A{3i)/B(31)/C(31)#Df31)»E{3H 

c #SI(162)#PHIR(162)/XC( lb 2 )>YC ( 162) #FM(162) # ARCL( 162 )• / i 

4 'rP4?3J5,^rP5{3Jj^‘^^^’'^‘^‘-°‘^^ 2)/ARC0L0(162),0EL0L0(162) ‘ 
COMMON /y PI»TP/FAO/fcM/ALP>ftN/PCH,XP/TC»CHO/OPHI/CL»RCL/YR 

»TK!jK*T?^TTY^^Jnn'^T-^P^' ^SE F/ SE PM/ TUE { 4 ) / M/ N/ MM/NN/NSP 
iNPTS# LL> I> LSEH^ «^^N£W#EPS1^ NC tS> XLEN# SCALOi 

/ INC/OFAC/GAM/KOES/PLTSZ/CFL/OPU 
Real mach/machn/new/Machs 

DIMENSION Hyiy), SEPpa62),CPP( 162)/ THETAP( 162 J/DELP( 162 J 

“iJ!fJ^^^^‘l»»DELS{l)/XX(l),Yy{l)/HACHN(l) 

E i(H(l)/FP(l/6))/(THETA(l)/FP(l/6)) 

Z , /CXX(1)>FP(1#3 ))^{YY(l)>FP(l#b))»(DELS(])#FP(l*10)l 

2/ (ANGNEw(l )/FP(l/24))/{SEPR(l)/FP(l/14n/(CPX{l)/PHlR(in 

I <'''A'^HS(l),FP(l/26))/(DSDm)/FPa/30>) 

4 >(DELX(1}/FP(1/12)),(TD{1),FP(1/20)) rr i i/.u; , 

CP(Q) a C5*((C4/(1,+C2*QA0 ) )**C7“1. ) 

USX(Q) a <C4-(1.+Q/C5***<1./C7) )/C6 
MACH(O) = SORT (0/ (C1-C2*Q ) ) 

Data ISW/0//CDF/0.//XPLT/.5//XFAC/100./ 

IF (NDES.GE.l) GO TO 5 
DO 1C J a i,NN 
PHKMM/J) a PHKl/ J)+DPH1 
PHI(MM+1,J) » Phi (2/ J)+OPHI 

y (ISW.EO.O .ANO.CSTAR.Ea.lOC. ) CALL GOPLOKhRN) 

COMPUTE AND STORF CP CRITICAL 
CPX{MM+1) a CP(1.) 

ISX + FSYMa3 IF FLOW HAS NUT BEEN CuMPUlEO 

I5X a {NCY+l)t(ITVP-3)*A6S{FSYM+iO. )+.2 
IF (ISX.NE.l) GO TO 30 
M4 a N3 
FSYM a 0. 

ALP a 0. 

XSEP a AMAXl (0./ XSEP-1. ) 
es a A(MM) 


20 


30 


DO 20 L 
X0L0( L) 
YOLD(U 
MaCHN(L) 
CPX(L) a 


IF 

GO 

DO 

L 

CS 


1/ MM 
‘ XC(L) 

> YC(L) 
a MACH(A(L ) ) 
CP(MACHN(L) ) 


( (ABS(YC(MM)-YC(l) ) 

TO 110 
40 L a 2/M 

« (PHKL + l, 1)-PHI(L-1/1) )*DELTH-SI(L) 

- {U*U)/FP(L/D 


.LE.I.E-5 ) .AND. (lAdS(NRN) .GT.999) » GO TU 50 


MACHN(L) a MACH(CS) 
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AO CPX(L) » CPChACHNU)) 

HACHN(MH) » .5*(HACHN(2)+MACHN(M) ) 

KACHN(l) » MACHN(MH) 

CPX(l) « CP(HAChN(D) 

CPX(HM) * CPX(l) 

QS«OSX(CPX(MH)) 

IF( (INOCO.EQ, 1) .AND. (FSYrt.EQ.7) ) GO TO 50 
IF (FSYH.EQ.6.) GO TO oO 
IF t (FSYM.LE«5») .OR* (ITYP.LE»2) ) GO TO 50 
C ADVANCE PLOTTER PAPER TO THE NEXT BLANK PAGE 

PL^T(12.0*FLUAT( INT{(20.2fXPLn/12. ) )/0.>-3) 

XP LT * #5 

50 CALL GETCP(CDF) 

IF(INOCO.EQ.l) ISW»1 

if(inoco.eq.i) return 

CALL GOPRIN < HP# THETAP# SEPP# CPP> DELP/ XTRANS ) 

IF (ISX.EQ.l) CALL EXIT 

ISW « 1 

RETURN 

60 DO 70 L « 1#MM 
70 CPP(L) = CPX(L) 

IF ( ( ISW.E&.O) .UP, ( FSYM.NE .£). ) ) GO TO 90 
C FIND THE BASE PRESSURE 
OELBP = 10. 

CPO « CP(hACHN(lXX-I n 
DO 80 L » IXX/M 
CPN » CP(KACHN(U) 

DEL3P a ANINl (OtLBP#CPN-CPU) 

80 CPO « CP^5 

BCP = BCP+RBCP*0EL8P 
90 ISW a 1 

PCH = ABS(PCH) 

IF (LSEP.GF.KM) CO tq j^xq 
C MODIFY THE MACH DISTRIBUTION 

CPO a CP(HACHN(LScP) ) 

SEPX a XC(LSEP) <{•> 


100 

110 


CM AX a OH IN 

DARC a TP/FL0AT(NPTS-1) 

DO 115 L a 1#NPTS 

115 HID a FLOAT! L-IH-DARC 
H(NPTS) a TP 

DO 116 L a i,n 

116 YY(L) a FL0AT(L-1)*DT 
YY{MM) a TP 

CALL SPLIF (MM»YY» ARCL»DSCT#C0#TD,3,0.,3,0.) 
CALL INTPL (NPTS»H#S# YY# AKCL>DS0r#C0>TC J 
S(NPTS) a ARCL(MM) 


SL a (BCP-CPO)/(XC(MM)-SEPX) 

DO 100 L a LSEP#MM 

CPP(L) a CPO+SL*(XC(L)-SEPX) 

MACriN(L) a MACH{QSX(C?P(L) n 

KUMIN a 1 

KOMAX a 1 

OMIN a MACriN(l) 
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CALL SPLIF IMM,APCL»MACHN»OSCT#LO»TD»3#0.,3»0.) 

CALL INTPL(NPTS/S »rtACHi/ARCL»MAcHN#OSOT»CO#TO) 

CALL SPLIf (MK>ARCL»XC»DS0T»C0> T0»3»0.»3»0.) 

CALL INTPL (NPTS#S»XX,ARCL>XC>OSOT#CC)»TO) 

00 120 L ■ ItNPTS 
IF (HACHS(L) .GT.QMAX) KQHAX « L 
IF (HACHS(L).LT.OMIN) KOHIN - L 
QhIN • A.^INKMACHSCD^OMIN) 

UMAX = AMAX1(MACHS{L)» JHAX) 

SEPR(L) » 0. 

H(L) - 0. 

OELS(L) » 0. 

120 THETA(L) » 0. 

IF (PCH.LT.O.l GO TO lAO 

KOHAX « KgMIN + lNOExR(PCH»XX{KOhlN*-l),NPTS-KQmH) 

IF (K0MAX.GE.NP1S) CALL A30RT 
lAO CALL NASHhC (KUF!AX,NPTS) 

XTRANS » PCH 

IF (PCH.LT.O) XTkANS * XX(KgMAX) 

KCiBOl « INDEXR(XTRANS>XX,KOMIN) 

IF {K0B0T.LF..1) CALL AJORT 
CALL NASHKC (KCbOT#l) 

FAC»S(A)/(S(A)-S(2)) 

THETAi 1)»FAC*THETA(2)+(1.-FAC)*THETA(A) 

H(1)»FAC*H(2) + (1.--FAC)*H(A) 

DELS(1)-H(1)*TH(TA(1) 

IFCXSEP.GE.O.) GC TO lAl 

FAC=(S(NPTS-3)-S(NPTS) )/(S(NPTS-3)-S(NPTS-l) ) 
THETA(NPTS)»PAC*THETA(NPTS-1) + U.-FAC )*THElA(NPTS-3) 

H(NPTS) »FAC*H(NPTS-1 )+(l.-FAC)*H(NPTS-3) 

OELS(NPTS) =H(NPTS)*ThErA(NPTS) 
lAl CONTINUE 

C COMPL-IE THE SKIN FRICTION DRAG 

0 = SOPT(OS) 

RI « (C1-C2*CS)/(C1-C2) 

htJT = (H(NPTS)+1. )<=(!. -C2*0S/C1)-1. 

HBB * (H(1)-H.)*(1.-C2*0S/C1)-1. 

CDF « 2.*THETA(NPTS)*0**(.b*( HBT ♦5.))«RT**3 
CDF « C0F+2.*THETA(1 .5*( HbB+b . ) ) *RT**3 
IF (ISX.tC.l) GO TO 200 

C HAKE OISPLACEHENT HuNJTJnE INCREASING ON THE UPPER SURFACE 

DO 170 L » KOMAX»NPrS 

IF (DELS(L*1) .LT.CELS(L)) DELS(L+1) - D£LS(L) 

170 CONTINUE 

C LUWEk SURFACE - FIND wHEKE DELS STARTS OECkEASlNo 

C TkEAT THE LOwFR SURFACE LIKE THE UPPER SURFACL IF XSeP.LT.O 

XPC » .60 

IF ( XSEP.LT.O.) XPC « 2. 

J = K 0 3 0 T 
180 J = J-1 

IF (LFLS(J-l) .LT.CCLS( J) ) cO TO 186 
U ( J .GE .2 ) CQ TO IdO 
GO TO 200 

185 IF (XX(J) .GT.XPC) GO TO IOC 
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1 


190 


200 


MUST STAY HONOTONb DECREASING 


Dfc’LS(J-l) « OELS(J) 


210 

220 


221 


230 


0ELS{J-1> « DELS(J) 

GO TO 180 
DISPLACEMENT 
J « J-1 

IF (DELS' J-1) .GT.DFLSUH 
IF (J.GT.2) GO TO 190 
SMOOTH DELS IS TIMES 
IF ( IS.L6.0) -GO TO *220 

00 210 I « 1>IS 
OLD > OELS(l) 

DO 210 L » 3»NPTS 
NEW » DELSIL-l) 

OELS(L-l) * .2b*(CLD+N£W*NtW+0ELS{L) ) 

OLD « NEW 
XPLT ■ XPLT+»5 

FAC«(S(NPTS-1)-S(NPTS))/(S(NPTS-1)-S( NPTi-2 ) ) 
DELS(NPTS)»FAC+0ELS(NPTS-2)+{1 .-FaC)*0ELS(NPTS-1) 

IFIXSEP.GE.O. ) GO TO 221 
FAC=(S(2)-S(l))/(S(2)-S(3) ) 

DELSm»FAC*0tLS(3) (1. -FAC) ♦Of LSI 2) 

CONTINUE 

IF (iSX.fO.l) GO TO 260 

Oh • (HIKQMAX+1 )-h(KCd.)T-l ) ) / FLOAT ( 2+KGMAX-K.CMlN) 

FAC ■ AKCCLO(NT)/S(NPTS) , 

IF ( XPLT* L T « 1» 2) CALL $ YMBOL ( • 35 >8 . • 1 55rii)IS*^LACEMENI 

1 AT EACH BOUNOAkY LA^ER IT ER AT lONf 270 .# ?5 ) 

CALL PLOT (XPLT + XFAC*0cLS<l)»lO.5>3 ) 

DU 230 L « 1»NPTS 

CALL PLDT(XPLT+XFAC*DElS(L)#10.5-YFAC*S(L).2) 

IF ( C L *GE .KUBOT ) . AND . ( L »LE .KQMA X ) ) H(L) = ri(L-l)+OH 
YY(L) = S(L)*FAC 

YY(NPTS) » ARCOLO (NT ) , 

OELX WILL 3E BOUNORY LAYER OISPLACLMENT AT NT POINTS 
CALL SPLIF(NPTS»YY#UELS/030T#CC#TU>3>0.#3»0. ) 

CALL INTPL(NT»AkC0L0»0ELX>YY/0fcLS#0SDT»L3»Tt‘) 

THE FOLLOWING ARE BEING COMPUTED FUR FLTuPE PRINT OUT 

CALL SPLIF(NPTS#S»0ELS#0S0T>CD>T0#3»0*»3fO.) 

INTPL(HM» AkCL»DELP»S»OELS/DSOT»CO#TO) 

SPLIF (NPTS#S>H#0S0I»C0»TD>3/O.»3>0.) 
INTPL(MM»ARCL»HP»S»H»DSL'T»CO»TD) 
SPLIF(NPTS»S*THETA»US0T»CO»TD>3»O.»i#O.) 

INTPL (HM»ARCL» THETAP»Sf THET A» OS DT / C0» TO ) 

SPLIF(NPTS#S» SEPR»DSDT»C0»Tl)#3#C.»3»0.) 

INTPL(MH»ARCL>SEPP#S#SEPR»DSDT»CO»TO) aitvTc 

THE SLOPES FUR THE O'^TEk AlKFOlL AT CUkRESPONtlNG POINTS 


frilCKNESS 


i 


CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

GET 


DU 2 AO L » 1»MM 

DDEL * ROLL*(DELP(L)-)So.‘*..U ) 

OELP(L) • DDEL 
DSUMIL) = DSUM;L)*D0£L 
2A0 SIL) * FAC*ARCL(L) 

S(MM) » ARCOLO(NT) 

CALL SPLlF(HH»S/FM»0S0r>C0»TD»3>0.»3#0.) 

CALL INTPL ( NT fARCOLD» ANGNEw> S»FM»DSDT/CO»TD) 
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bELHAX < 0. 

DO 2i>0 L « 1»NT 

DDEL « DELX(L)-CELOLO(L) 

DELHAX a AHAX1(DELMAX> A8S(00EL) ) 

DY » DELOL£)(L ) + RDbL*DOLL 
AN6 ■ .5*( ANG0L0(LH-AN6N£W( D) 

XX(L)«XOLD(U 
YY(L)»YOLO{L)+OY/COS( ANG) 

250 OELOLO(L) » OY 
ISS « IS 
IS a -1 

IF (ITYP.EQ.99) CALL GJPRIN ( hP» THEY AP> SEPP» CPP» OELP» XTRANS ) 
CALL AIRFOL 
IS « ISS 
FSYM » 7. 

RETURN 

260 DO 270 L = 1»HM 

ARCOLD(L) = ARCL(L) 

CPP(L) a CPX(U 
270 ANGOLD(L) « FM(L> 

CALL SPLIF(NPTS,S/OELS/OSOT>CU>TO#3>0.,3#0.) 

CALL INTPLtMrf^ ARCL» DSUM» S> DELS> DSDT» CU»T0) 

CALL SPLIF(NPTS»S»SEPR»DSOT>CC»TD> 3»0.»3»0. ) 

CALL INTPL (KM#ARCL»SEPP#S»SEPR#DSDT>CC»TD) 

CALL SPLIF (NPTS»S# ThETA#DSDTiCC>TD» 3»C.#3»0. ) 

CALL INTPL {NH»ARCL>THETAP>S>ThETA>OSDT»CO,TO) 

CALL GOPRIN (HP»THtTAP»SEPP>CPP»DELP>XTRANS) 

NT a MM 

CALL GETCP(CDF) 

IF (JK.LE.-l ) CALL PLOT (0./0./999) 

CAIL EXIT 
END 


SUBROUTINE GOPRIN (H» THETA# SEP# C PP> DEL/ XTR ) 

REAL KACHN 

COMMON PHl(lt2#31)#FP(162»31)#A{31)/B(3i)>C(31)#0(31)#E( 31) 

1 >RP(31)»RPP(31)»R{31)»RS(31)#Kl(31)/AA(i62)»BB(lt2)#C0(lb2) 

2 #SI(162)»PHlR(lo2)»XC(lo2)#YC(162)» FM(162)» AkCL(162)»DS0M(lc2 ) 

3 # ANGCLD( 162)# XOLC( 162)# YOLO (162)# ARC OLD (162)# OELDLD( 162) 

A #R?A(31)#RP5(31) 

COMMON /A/ PI#TP#PAO#EM#ALP#PN#PCH#XP#TC»ChU#L'PHI#CL»RCL.YK 

1 #XA#YA# TE#DT#UP#DELTH#DELR#KA#DCN#0SN#KAA#EPSIL#QCR1T#C1#C2 

2 »C4#C5#C6f C7#&ET#BFTA#FSYM#XS£P#SEPM#TTLE(-^)#M#N#MM#NN#NSP 

3 #IK# JK# 1Z#1TYP#H0DE# IS#NFC#NCY#NRN#NG# lOlM# N2# N3#NA#NT# IXX 
A #NPTS#LL# I#LSEP#MA#NEW#CPS1#NDLS#XLEN# SCALOl 

5 # SCALQ0#N6#GAMMA#NQPT#CSTAR# Kf M# DEP# OINF# TSTEP# XUUT 

6 # INC»QFAC.GAM#KCES#PLrSZ#(jPL»UPU 

DIMENSION DSDT(1)#KP'>(1)#FPPP(1 )#ri(l )#SEP(1)#THETA(1)#CPP(1) 

1 »MACHN(1)»CP{1)#DEL(1)#BL(A) 

EQUIVALENCE {FPP(l)#cJ(l))#(FPPP(l)#SI(l))#(DS0T(l)#FP(l#31)) 
equivalence (MACHN(l)# A(1))#(CP(1)»PHIR(1)) 
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DATA l0N/l0FF#Z#SEPHAX/l>0/0.# ,004/ 

SN ■ -2./ARCL(MH) 

QHIN«NACHN(1) 

DO 10 L « 1>H 
QHINoANiNKHACHNUUOMIN) 

10 ARCLCL) ■ AC0S(1. + SN*ARCL(U) 

ARCKMMI « PI 

CALL SPLIF(HM#APCL»FH»OSOT#FPP, FPPP> 1/0.»1»0, ) 

OSOT(l) • FPP(l)*l,E-5 
DSDT(HH) ■ -FPP(HM)*l.t-;» 

DO 20 L » 1>HM 
FPP(U a RA0*FM(L)-180. 

20 FPPP(L)» SN<'0S0T(L»/A«AXU1.E-!»,SIN(ARCL{L) ) ) 

IF (FSYH.GT.S.) GO TO 120 
IF (FSYM.EQ.O.) GO TO &0 
WRITE (N4,310) 

2tf IF (FSYM.EQ.O) WRITE (.^^4>320) TILE 
IF (XP.EQ.O.) WRITE (N4»3o0) lOFF 
IF (XP,NE«0.) WFITE(N4>660) IDFF 
IF (XP*E0,0,) GC TO OOO 
REWIND H4 
READ 

READ (M4#5i»5J 

555 FORMAT (IHl) 

6C0 DO 30 L » 1»MM 

IF (XP.EQ.O.) GO TO 610 

IF (HOD(L+l#55).£Q.O) WRITE (N4/660) ION 
PEAD(M4,556) OUM.OUM 

556 F0RHAT(2F10.4) 

WRITE(N4,670) L» XC ( L )> YC ( L ) # F PP ( L ) , F P PP ( U » MACHN ( L ) / CP ( L ) > OUM 
GO TO 30 

610 IF (MOO(L+l»55).EG.O) WRITE (N^>360) ION 

WRITE(N4>260) L» XC ( U > YC ( L FPP ( L ) » F PPP (L ) .MACHN ( U »CP ( L ) 

30 CONTINUE 

670 FORMAT ( 114# 2F9 . 5» 2FS . 2» 3F9 ,4 ) 

C RESTORE QUANTITIES TO VALUES THEY HAD UPON ENTERING THIS ROUTINE 

40 DO 50 L » 1#MM 

ARCL(L) » (C0S(ARCL(L))-1. »/SN 
50 FP(L#NN) « 1. 

CALL COSI 
RETURN 

60 RNX « .1*AINT (PN+l.E-5) 

IF ( ( ABS(YC(MM)-YC( 1) ) ,LE .l.E-5 ) . AND. ( 1A3S(NRN) .GT.999) ) GO TO 25 

WRITE (N4»390) TTLE»RNX 

WRITE (N4#330) lOFF 

IF ( JK.GE.O ) GO TO 80 

CALL PLOT (2,#0.,-3) 

ENCODE (30>37C#TTLE) EM#CL#TC 
CALL SYMBOL ( 1 .2# ,7# , 14# TTLE# 0.# 30 ) 

ENCODE (20#380»TTLE) RNX 

CALL SYMBOL ( 1 .5, 1.0# , 14# TTLE# 0, # 20 ) 

CALL PLOT{PLTS/*XCll)#t». + PLTSZ>!‘YC(l)# 3) 

DO 70 L •= 2# MM 

70 CALL PL0T(PLTS2*XC(U#3. + PLTSZ‘i‘YC(L),2) 
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IPEN « 3 

00 DO 100 L ■ 1#MH 

XS » X0L0(L) + 0SUM(U*SIN(ANG0L0(L) ) 

yr ri* 1 ^ ^ I- * ^ 

*C(L) • XS 

YC(L) » YS 

IPEN‘^!i*“*‘^’ '’‘•'^^“’‘•T52*XS,5.*PLTSZ*YS,IPEN) 

IF JM00(L+3#55) .EQ.O) ^RlTt (NA,330) ICN 
IF (XOLDCD.GT.XTP) GO TO 90 
TRANS » IH 

Tc TRA.-JS * lOHSTAGNATION 

IF GO TO 85 

*“'-““’'''“‘-“"-''fP>'‘Ll.fPPP(Ll,CPP(U,THErAlL) 

100 CONTINUE 

IF (XP.EQ.O.) NRN « -IABS(NRN) 

XP « -ABSIXPJ 
RETURN 

120 WRITE (N<f,310) 

WRITE (N4#300) lOFF 
I ■ 1 

YSEP « ASS(XSEP) 

IF (XSEP.6T.0.) YSEP » 2. 

DO 150 L » 1,«H 

\l .EO.O) WkITE (N4/300) ION 

IF (XC(L).GT.XTR) GO TO 130 
TRANS « IH 

TP f TvriMl'i TRANS » lOHSTAGNATION 

IF^(UC(L+1).6T.XTR).0k.(XC(L-1).GT.xTR)) go to 125 

YSEP a ABS(XSFP) 

GO TO 150 
130 bL(l) a IH 
BUZ) a IH 
BL(3) = IH 
EL(4)a IH 

IF (L.EO.LSEP) BL(1) * 2HLS 

BL(4) = 2HLP 

WRITE (N^,280) faL ^ L# XC ( L ) > YC ( L) / F PP ( L) , F PPP ( i ) , HACHN < l 1 
GO TO 40 

260 F0RHAT( H4,2F9.5»2F8.2»2F9.4) 

^'^^IEQ^Z^IsT' Pe. 2*F8. 4> 2F9. 4, F9. 5, F9. 5, FV.5>F7. 2, 
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290 FORMAT t H6>2F9 »5# F9.2» F8 ♦ F8 • A# 2F9»4#2F9» 5» 3X» A10#7X# I i 1 

310 format ( IH l# I jX# 40HL0WER SURFACE TAIL TC UPPER SURFACE TAIL ) 

I!a 17X26HLIST1N3 OF CODROINATES F0R,2X,4A4) 

330 ( II /llXlhX»8X» IHY^ 6X# 3HANG# 4 X j ShKAPPA# 6X» 2hC P» SX# 

1 5HThETA,SX,3HStP,6X,2HXS,7X,2HYS/) 

|P1^»!’»F9.5»Fo.2,F3.c,F9.h>4X,A10>4X>2F9,5J 
350 FORMAT ( F14* 5> F9.5> F8 .2> F6 .2> F9 , 4# 4F9 , 5 ) 

? 7 n TnoMfl J 6X,3HASG,4X5riKAPPA4X4HHAv.hcX2HCP/. 

370 FORMAT ( 2FlH», F4 ♦ j, 4X> 3 mCL => F5 . 3#4X» 4HT /C «» F4 *3 ) 

380 FORMAT (4H RN*>F4.1>9H MILLION ) 

390 FORMATdHl/ 9X26HLISTING OF COOKOINATES F0R#2X#4A4 »4X#3riKN»/ 

1 F4.1»8H MILLION ) 

I I1^12X1HL»6X»1HX» 3X> 1FIY» 6X# 3HANG>4X5HKAPPA4X4HMACMcX2hCP» 

X oX^4HDATA/) 

END 


SUBROUTINE NASHKC (R1/K2) 

COMPUTE TFIE BOUNDRY LAYER FP.Uh POINT K1 TO R2 
K3 WILL BE THE SEPARATION POINT 

COMMON PHI(162»31)»FP(162f31)^A(31>»o(31)»C{31)>0(31)»E(3i) 

f ^^.^![i^^i^^^^'^*^^^^*^^*^^2)^YC{lU2)»FM{lo2)»AKCL(ic2J#D5UM(162) 

3 ^AN60LD(l62)>X0LC(l62)»Y0L0(iU2)»Ar<C0L0(lo2)/i;ELCLC{lb2} 

4 /RP4(31)»RP5{31) 

COMMON /A/ PI#TP»RA0#EK* LP;>RN, PCH# XP#TC»CriO*DPhI#CL#RCL# Y'< 

1 »XA^YA»TE,DT,Dr»0ELTH,0ELR,RA,DCN,DSN^RA4,E?SlL»CCklT,ti,U2 

I (4),M,N,MM,Nt.,NSP 

1 * ioTt^i }^r NCY> NRN. NG» 101 M» N2/ N3# N4# N1 ^ 1 XX 

4 »NPTS»LL#I»LSEP»M4fNcv»»EPSl^FDES^XLEN>SUAL0i 

5 #SCAL00#N6#t>AMMA>NQPT#CSTAR>R£M^0EP#QINF*rsrEP»X0UT 

6 >inc»qfac>gak,koes.pltsz,qpl>cpu 

DIMENSION MACHS { 1 ),H{ 1 ), theta m,SEPR Cl J,S{l),DELi(l),XX(l) 

EOUIVALENCEC MACHS (])>F?(1j 2B))/(H(1)/FPCI/6))»(T He TAClJ^FPClic)) 

REAL MH»MhSU,NU>MACHS 

DATA TR,»TH0,TEi,TE2>S£PMAX,PIMI,N,PlMAX / . 34 24, 320 . , 5 . £-3, & . L--b . 
1 .004,-1. 5, 1.E4/ 

GAMI = ,5/C2 

CSIINF » C4 

INC * ISIGN(1,K2-K1 ) 

YSEP » ABSexSEP) 

IF ( CXSEP.GI.O.) .ANU.CINC.LT.O) ) YSEP = 1. 

SEPMAX s SEPM 
GE = 6.5 
L « K1 

DS » AaS(S(L)-S(L-INO) 

' LP = L+INC 

MH « .5*CMACHS(L)+MAChS(LP)) 
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HhSJ - iiLPRODUClBILlTY OF THE 

CSIH » l»+C2*MHSO ORIGINAL PAGE IS POOR 

OSOLD « DS 

OS » A6S(S(LP)-S(U) 

PHOH » T^^GAIil 

NO » T*(l.+T8)/(RHCh«{T*TR)) 

RTH « RN*MH/( EP.+NU) 

IF (L.NE.Kl) 60 TO 3C 
THETAHs RTHO/PTh 
THT a THLTAH 

30 FC * I *0+ •066*tthSC-. C03*Mri*MHS0 
_ FK a I •“•134*MHS0+.027*MHSQ*NH 
- DU AT MOST 200 ITERATIONS 

DL i<*0 J a 1»^99 

TAU a RTAU+RTAU 
HB a l./(l,-GE*PTAU) 

HH a (He + l.) + (l.*.i73*N,Hso)-l. 

SfcP = -ThETAH*COOS 
IF (S6P.LT.SCPMAX) GL TO 50 

t,n SEP a SEPHAX 

^0 PIE a HH+SEP/TAU 

PIE a AMAXKPIPIN, amInKPLMAX, PiE) ) 

G a t.l*SJKT( PIt*1.81)-l, / 

T2 a ABS(6-GE)/6E 
GE a G 
0T2 a OTl 

CTl a (HH+2.-KrtS0)»iEP+TAU 
IF (J.EO.l) GO IC lie 
T1 a ABS({DT1-DT2)/DT1) 

no TO 130 

no THETAH a THT<-.5*DT1*DS 

140 CONTINUE 

130 THETA(LP) a THT+DT1*US 
SEP a -THETAH*DCDS 
THETAH a TriETA(LP) 

THT a THETA(LP) 

lEPR^‘L J) \ »'‘‘^SfSEP*USULG ) / (OS+DSOLC) 

niLP/ ~ HH 

CcLS(L) a H(L )*THETA(l) 

L a LP 

IF (L.NE.K2) GO TO 10 

oc.c " 2.*SEPP(K2)-SEPMR2-INC) 

0ELS(K2) a H(K2)*THETA(K2) 

H(K1) = 0. 

SEPR(KI) = 0. 

CALL NASHLSC<2) 

DELS(K2-INC) a H ( K2-I NL ) *Tht ( K 2-1 NC ) 
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0E^LS(K2) * HtK2)«THETA(K2) 

RETURN 

END 


SUBROUTINE TRIOK A»B»C» RHS»OUT#N> lOIN) 
COMPLEX RhS/OUT 
DIMENSION A(l),e(l),C(l) 

DIMENSION GA{ 35 )> RHS ( ID IM ) / OUT( 31> ) 
REC=»1,/B(1) 

GA(1)»REC*C(1) 

OUT(l)=iREC*RHS(l) 

DO 10 J-2/N 

REC«=1./(3(J)-A(J)*6A(J-1) ) 

6A(J)»REC*C( J) 

10 0UT(JJ=REC*(RHS(J)-A{J)»0uT(J-1) ) 

DO 20 JJ>2>N 
J=N-JJ+1 

20 CUT(J)=0UT(J)-GA(J)*0UT{J+1) 

RETURN 

END 


1 


SUBROUTINE SOLVl 

COMMON PHI(162>31)>FP(lo2>31)>A(31)»B{31)>C(31)»0(31)»E(3i ) 

1 ^RP(31),RPP(31),R{31),RS(31),PI(31),AA(162J,6B(ie2»,CG(162) 

? *SI{162),PHIR(162),XC(lo2),YC{162)^FM(162),ARCL(162j>0SUM(l62) 
3 » ANG0L0( 162)# XCLO( 162) >YQL 0(162 )> ARC OLD (162)# DElOLOU 62) 

A #RPA(31)#RP5(31) 

COMMON /A/ PI>TP»PAO#EM#ALP#RN»PCH#XP#TC#Chu#oPHI#CL#RCL#YR 
I # XA# YA# TE#DT#DR#DELTH#OElk#RA# DCN# DSN#KAA#hPSlL#CCRIT#Cl#C2 
f CA#C5#C6#C7#8ET#BbTA#FSYH#XSEP#StPM#TTLE(A)#M#N#MM#NN#NSP 

# IK# JK# 12# ITYP#MODt# IS#NFC#NCY#NRN# NG# I31M#N2#N3# N/f#NT# IXX 

»NPTS#LL#I»LSCP»M^#NhW#EPSl*NDES#XLEN»SCALOI 

# SC ALQ0#N6#GAHMA# NQP r#CSTAR#FbH#iJEP#QlNF#TSTEP#XQUT 
# lNC#QFAC#6AH#KDbS»PLTS2# OPL#CPU 


complex FF#F1»G6 

DIMENSION CX( 162)#SX(i62)#FF(lo2)#GG(162)#ri(31) 
COMMON /SOLI/ 0(162#31) 

IF(NbW.NE.l) GO TO 30 
DO 1 I«1#M 
CX(I)=COS( (I-1)*DT) 

SX(I)=SIN((I-1)*0T) 


CONTINUE 


30 


NEW=0, 

MMP»MH+1 

CONTINUE 

Ma=M/2 

MA1=MA+1 

DU 2 J»1,N#2 



DO*’? 

IH-M-I+3 

o( r» j) « real(ff( I ) j 

<2{I# J + l)=PcAL(G6(I)) 
Q(IM»J)=-AIflAG(FF(l)) 

0( ^M» J+1)=-«IMAG(6G( 1 1 ) 
7 CONriNUfc 
2 CCNi'InUE 
HRs , ‘ +DR 


DO 3 J=l,n 
D( J ) =2 .*rs (J ) 
T»RA*RA*R( j) 

3 C(J)=T*{R( J)+hR) 

ca)«o(i) 

DC <* 

IH«M-I+3 
DO 5 Jsi,^ 


a«J>=-D(J)-24*(1.-CX(I)) 

CALL TRI01(B> A»C/FF#F1#N>16: 
DO 8 J*lfH 


0(I»J)sREAL(Fl»J)) 

0( IH» J) =AIf1AG(Fl( J) ) 
8 CONTINUE 
<► CONTINUE 


Du ■) Jsi,M, 2 
DU 10 Isi,fiAl 
IN = M-I + 3 


) 


10 

9 

12 


FF I .CMPLX(Q(I,j,,-o(i„,j,) 
?’:CNPLX(C(I,jn),-J(lN,j + i,, 

FF lN)sCMFLX(Q{I,j),c(i:i,j)) 

cSNriJiE 
CD 12 Jsl,N 
0(MH> J)=ca, J) 

0 ( H^P> J ) sQ (2» J ) 


RfcTUFN 

END 


SUBROUTINE SWCfPl 

lC0NHDNiPriI(l62,3n,FP(162,31),A^ 

2 » SI( 162 ), PHIR(162),XCf ?i v-f 1 • ^ ^‘-■^’*D0(162 ) 

3 /ANGDLDC 62),X0L0(l62),V0lD(ls3)‘^ARf t/om 

^ >KP^(31 ),RP6{31) ‘•^‘^‘’‘^*^ARL0L0(l62)A0ELuLD{lb2) 
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3 » IK> JK»IZ#IIYP#M00E>Ii#NFC>NCYiNRN»NGiI0If'.».N2iiN3#N<.»NT»lXA 
#NPTS>LL/I»LStP»K<»*Nc«*>tPSl> KJtS»XLfcN»SCALCii 

5 *SCALQ0»Ne),GAMMA»NuPT»C5TAk/KfcM> 0bP>01Nf»l$TbP#XDUT 

6 » INC»0FAC»GAM,KC£S,PLT5>Z»QPL»aPU 
COHMUN /SbLl/ C(162»31) 

DATA 0/5022*0.C/ 

YP»0. 

NSP»0 

DO 10 J>1»NN 
PHI(MM»J)«PHI(1,J)+0PH1 
PHI{PH+i» J)«PHI(2» J»*bPHI 
10 CONTINUE 
TE»-2 

DO 30 1»LL#«M 
CALL HURNANl 
DO 100 J=1»N 

100 CONTINUE 
30 CONTINUE 
TE-2 
1»LL 

80 I»I-1 

CALL NURMANl 
DO 60 J»1>N 
0(I> J)>0(J) 

60 CONTINUE 

1E( I.GT.2) GU TO 60 
DO 61 J=1,N 

61 C(l> J)«£I(HH* J) 

210 F0RMAT(t)(2I^»flt,en 
CALL SOLVl 

200 F0RHAT(6( H»L16.8 n 
DO 110 1*1»M 
DC 110 J«1»N 

lie PHI(l.J)»PHni>J)+0(l,J) 

DO 111 J=1»N 

Hi PHI (PM> J ) *PH1 (1» J )+OPril 
1F(RCL.EO.O.) OL TO 90 

YA«RCL*((PriI(M»l)-(PHl(2/l)<-rPHI) ) ^ Ofc LT H + S 1 ( 1 )) 

IKHQDE.tO.l ) GD TU 90 
IF (NDES.GE.O) CU TO 9l 
ALP=ALP-.6*YA 
GO TD A2 

A1 B6(l) « BB(l)-.t*YA 
9Z CALL COSl 
GO TO 95 

90 YA*TP*YA/ (l. + BET ) 

DPHI-DPHI+YA 
95 DO 97 L»1»N 
97 PH1(L»NN)=0PH1*PHIR(U 
IMMCOE .Ew.O) rLFoftN 
DU 120 J«1,N 
DO 120 L=1»M 

120 PHI { L / J ) =PH1( L* J ) ♦YA*PH1P( L) 



RETURN 

END 


SUBROUTINE HURNANl 

CC'IMON Phl(162/31)»EP(lto2/31J/A{3D# s(31)#C(3i),0(31J»L(31) 

!• #RP(31)»KPP(31)#R(31).R^(31)#fvU31)/ AA(162),n3( ic2)>CO(162) 

2 jSI(162)#PHI«(lt>2)»XC(io2)»YC(162)»r".(i02)» AKCL(lfc2)*CjuM(l52J 

3 » ANGCLD( 162) » XQLD( lo2 ) » YGLC( l62 ; . AkCQL3(lo2 ) » DE LULC 1 162 ) 

A ,aP^(31)/RP!?(31) 

CuH.iL-N /A/ P1,TF,RAD> t.'l,ALP^RN,PCh» X) »TC»ChO»i)PHl»CL»RCL» YR 

1 ,XA, YA, Ic/OT/DPf DELlH,DbLK/PA>CCN/ JSN/RA4» =PSU/U(,k1I#C1#C2 

2 »C<»»C5»C6»C7>PbT»bEl A»FSYM,XSLP»StP.'1»TTLb(«.),N»N,N't»NN»NSP 

3 » IK» JK. IZ, ITYP,miDtf Ii»i\FC>NCY»NKN> NG. I01N.N2»N3»N4,NT, IXX 

4 /NPT S»LL# I»LSLP#N4,NEw.. tPSl/ XLbN. SCALUl 

t> , SCALC0»N6»UAKr.A#NCFWCSTAK> ktK.» 0tP#01\F» ISTbP# XuUT 
6 t INC# 3 FAC»GAM»KuFS»t'LTiZfQPL^'JpU 
PHIU = PHI( l,2)-2.«f'R*CJ( I ) 

PHlYP»PriI(I,2)-PHl(I#l) 

PHlYY=PHIYP+PhlL-PHI (1,1) 

PH1XX=PHI( I*1,1)+PHI(I-1,1 )-PHi ( 1, i)-P^l( 1,1) 

PhIX««PH 1( !♦!, D-PHK i-1,1 ) 

PhIXP»PHI ( 1*1,2 )-Ph 1(1-1, 2) 

IKI .NF.MM) GO TC lu 

0( 1) =Cl*(PrilXX*Ki (1)*PR1 tY*fvA4*C0( 1 ) > 

D(1)«-D(1)/C1 
GO TO 40 

10 U=PHlXM*OELTh-Sl (I ) 
bC»U/hP(l,l) 

OS = U*iJO 
J»1 

IF(3S»LE.QCRiT) GO TO 30 

D(1)=0. 

GO TC 40 
30 CONTINUE 
CS=Cl-C2voS 

e« = 6C*cS*(FP( I-l, 1)-FP( I*l,in 
X = RA44(CS+tS)*Clj(I) 

CNOS-CS-JS 

U(l)=CS«RS(l)*PhIYY + H(i>»DO*X+CMJS*PMl\K 
L( 1 ) =-0( 1) /CS 
40 CONTINUE 
OL 6C J=2,N 

PHIXX.PHI ( 1*1, J)4FH1 ( I-l, J )-PM ( 1, J )-PHI ( 1, J ) 

DU=PH IXP 

PHIXr=PHI( I*l,J*l)-PhI( I-l, J*l) 

PhIXY=PhIXP-PHlx^ 

PHI Xb *0U 
OU*DL*DtCH 
PHIYb *PHI YP 

PH1YP=PMI(I,J*1)-PH1(I, J) 

PHI YY^PHIYP-PHIYN 
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U*R(J)*OU-SI( I) 

RAV«R( J)*r<A*V 

SQ«1./FP(I,J) 

BOUaBO*U 

'JSa 8 CU*U 

UV«( BQU+BOU)*V 

VS«BO»V*V 

OS-US+VS 

IHQS.LE.QCRIT) go to 50 
0 ( J ) *0 • 

60 TO 60 
50 CS»C1-C2*0S 
CMVS«CS“VS 
CK'JS«CS-US 
UVla,5*BQU*RAV 
C(J)«RSfJ)*CKVS 

50 CONTINUE 
return 

END 


L - N/2 

iI\kI complex TRA^FORM 

if (NS, LI .0) GO TO 20 
00 10 J • l,N 

Gc'u.’.o 

JO jU'.i'’ '‘'’‘"S' TkXNYfotn 

DO 30 K « 1,1. 

>'XI‘'^XO(XLP(K) J+ktAL(BcT(K 

30 J a 


i 
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<"» o 


K-l+1 

«0 CALL FfO«(UN,AlA,BtT,CS,SV. REAU6E 

no; separate out the R^I^'^Lo TRAOTHRR, pbrts 

(NS.LT.O) 60 • 60 

EM».5 

00 50 K « 1,L 

6C 00 7C J . l,N 

70 ^,4"'!’ * REAL(ALP(J))/N 
70 G(.,J) . -A1MAG(ALP(J))/n 
RETURN *o/i/r< 

END 


10 


ru.^,iJN VLAYER(EM2 ,a,£ 
X»{ bP2-A) / ( B-A) 

VLAYER « 0. 

IE (X.Lb.O.) RETURN 
IF (X.GE.l,) 60 TO IG 
VLAYEK a 3 •♦X^A— 2 • 


RETURN 

VLAYER 

RETURN 

END 


1 . 



iUBROUTiNb NASHIS(K2) 

PREVI00s''5''vAiut” '"’'^■^’’”''‘2),StPP(K2-l),SEPR{K2) 0SU(, 

^ #RPA(31)>RP5(31) ” It-*? )# DELOLl/( 162 ) 

COMMUN /A/ PI> rF>RAD#t<1#ALP.KN.PfH yp rr •i.r. 

1 >XA,YA, TE, OT,CR,OCLTri, OP! 

2 >C^>C5,C6,C7,BLT,BlTA.FsoIx\;^‘1'Sm 

3 'I*^'JK/ I2#irYP>/lOGE>IS#NFC»NCY i^ciu ^ ^ ^ NN» NSP 

^ '•'^PISpLL^IpLSc-PpMARNt.'^tPsI ' 

>SCALQ0.N6,GA«hA,NGRTrCi?AR 
6 'INC,OFAC,GA,‘.,KOES,P^iLJp, oSn 
OlnENSlON 3(l),H(l);s^rn 

^ IVALE.NCE ‘S(l),FP(l,i6,,,(H(l,,FPa,o)),(SEPR(x),FP(l,14)) 
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Aiii+s+s ♦ cm 


10 


F(S» 

NLS « 5 
XNLS s FLOAT(NLS) 

LLS « NLS + 1 
XI » 0. 

X2 ■ 0» 

X3 « 0. 

XA » 0. 

DO 10 L=2,LLS 

CUM « S(K2-FL0AT(L*1KC)) 

XI « XI ♦ OUM 

X2 a X2 + 0UH*0Ui'1 

X3 a X3 + DUM+CIM+OUN 

XA « X^ + UUN*CUK*i}UN<'OUH 

2211 s X2-X1+X1/XNLS 

2321 « X3-X2»X1/XNLS 

2A22 = Xt-X2*X2/XNL> 

Y1 a 0. 

Y2 a 0, 

Y3 a 0. 

DO 20 L=2>LLS 


20 


DOM . 
DOOM 
Y1 a 
Y2 a 
Y3 a 

22 a 

23 a 
Aill 
Dill 

cm 


30 


' S(K2-FL0AT (L*lNCn 
a H(K2-FL0AT(L*INC) ) 

Y1 + DUMM 
Y2 ♦ DUMM+CUM 
Y3 ♦ OUMM^DUM+UUM 
Y2-Y1*X1/XNLS 
Y3-Y1+X2/XNLS 

321*2^21) 

' (22-Am+2321)/Z211 
a <Yl-Alll*X2-Bm*Xi)/XNLS 
H(K2-INC) a F(S(K2-lNO) 

H(K2) a F(S(K2)1 

Y1 a 0. 

Y2 a 0. 

Y3 a 0. 

DO 30 La2,LLS 

DUH a S(K2-FL0AT(L*1NC ) ) 

DUMM a SEPR(K2-FL0AT(L*INC)} 

Y1 a Yi + 

Y2 a Y2 ♦ DOMM*OUM 

Y3 + 0UMM*DUM<-0UM 
Y2-Y1*X1/XNLS 
Y3-Y1+X2/XNLS 

a 3*22 11-Z2 *2321 )/(ZA22<‘Z211-Z 32 1*Z321) 

a (22-Alll*Z321)/Z211 
' tyi“AlU*X2-3Ul*Xl)/XNLS 
SEPR(K2-INC) a KS(KZ-INC)) 

SEPR(K2) a F(S(K2)) 

PE TORN 
END 


Y3 

22 a 

23 a 

Alll 

8111 

cm 
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RtAL NtW/LEFT ^ ^ ‘ i J ' PPPP( 1) , XI (1 )j 

data TOL /1.E->) / 

C ^*^P»FPp!^FPPP^ARrTHE*FUNCTTr^ FOR L • 1 

NX »IABS(MX) ^NCTICN aND DERIVATIVES 

K m 2 
LI » I 
NEW « X(l) 

FN » F(l) 

PVAL . FI(1) 

If UBscfm-fiun.Gi.Tou go to 5 

XHl) . x(l) 

5 DO 

GO TO 0 

FN s F(K) 

IF (fP{K)»FP(K-l).6T.O.) GU TO 6 

SGN > P(K-I)-FVAL 

Jo ’ ‘^N-FVAL 

IMFP(J)*fP(J-1,.le.u.) go to 7 

GO J^T^;^;<'‘''»-f=VAL).LE.O.) GO TO 20 

^ ox‘\^ -2!iFp\‘S!i!/7JpJuli;!cJ‘‘': 

RIGHT » X(J-l)+[,x 1 ^'*^SIgN{R0l:T>FPP{J-1))) 

IF (LEFT. GT, RIGHT) 60 TO 10 
F2 « .b*FPP(j-i) 

F3 = FPPP(J-1)/6. 

IF (MX.GT.O) GO TO 11 

MX » L-1 

RETURN 

11 PRINT A99,l,FI(U 

A99 format ( * trouble AT *, 17, 3X, F 1 6 . 6 ) 

GO TO 100 

20 OLD = am AXl (X ( J-l ),|r,Ey+]-yi^ j 
F2 » .i>*FPP(j-i) 

F3 » FPPP(J-1)/6, 

START-OLD 
DO AO K - 1,10 


FI(l) 

TU ABS(MX) 

AT THE X POINTS 


1 )/ 6 . ) ) 
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ox ■ OLD-X(J-l) 

IF (ABS(FPOLO) .L£»TOL) 60 TO 60 

FN » F( J-1 )+0X*( FP( J-l ) *ux^{ F2+DX*F3) ) 

NEW » 0L0-(FN-FVAL )/FPJLO 

IF (NEW. LT. START) GO TJ 60 

NEW » AfTNl(NEW#X(J)) 

IF t ABS(NEW-OLO).LT.TJL) GO TO 90 
40 OLD > NEW 
CALL ABORT 
60 RIGHT ■ X(J) 

LEFT » OLD 

IF (SGN»(FN-FVAL) .GT.O.) GO TO 65 
RIGHT « LEFT 
left • XI(L-I) 

IF (L.EO.l) LEFT * Xvl) 

65 DO 70 K 8 1>50 

IF ( (RIGHT-LEFT). Lt. TOD GO TO 90 
NEW « .5*(LEFT+R16HT) 

OX ■ NEW-X(J-l) 

FN . F(J-l)+0X*(FP(J-l)t0X*(F2+0X*F3)) 
IF ( ( FN-FVaL) ♦SGN. L£ .0. ) 60 TO 80 
LEFT» NEW 
60 TO 70 
80 RIGHT ■ new 
70 CONTINUE 
90 XKL) » NEW 
100 K » J 
NX » NX 
RETURN 
END 


SUBROUTINE REAOOS 

PRESSURE DISTRIBUTION 

COMMON PHI(162#31)>FP(162#31)»A(31)#ti(31)>C(31)#0(31),E(31) 

2 »SI(162),PHIPa62)f XCll62),YC(l62),FM(i62),AKCL(162)»LSufmr^2) 
I 'rSJ“ 3 J[,^rp’J 5 JJ^‘^^^‘'^‘^‘-‘'‘^^ 2 ),ARCOLU( 162 ),uElOLL( 162 ) 

COMMON /A/ PI#TP,RAp,EH,ALP,kN,PCH,XP,TC,CHD,OPHl,CL,RCL,YK 
1 »XA^ YA» TE» OT. OR. C-ELTH»O£LR^hA.DCN.0SN^RA<»»EPSIL#QCHT»Cl»C2 

1:>>NFC#NCY,NRN,NG# IDIM,N2,N3.NA,NI,IxX 
>NPTS»LL#I»LSEP#M9.Ncrt>cPSl>NDES>XLEN^SCALQI 
#SCAL00.N6»GAHMA'NQPT>CS FAR/ kEM.OEP, 0INF> TSTEP> aOUT 
»INC»QFAC»GAH»KBES»PLTSZ»QPL.OPU 

•‘^‘^’'5F‘1>'9X(1),SX(1),ES(1),GP(1),GPP(1),GPPP(1) 

1 >D3DS(1)»PHT(1)>DPHOS(1).6PP{1)^UPPP(1),0(1) 
equivalence (0K1).FP(1,1) ),(SF(1),FP(1,3)),(0x{1),FP(1,5)) 


156 


?>F Poca glUAUTY 



500 


510 

20 


30 

bOO 


60 


ITATlnL'x $ 3 “ 2 j;"'‘-'^’‘''*""*'l-'*“rA<.CSUR-.i,.<OAI,. 1 A-l 

MODE > 0 
REWIND N6 

REAO{N6/500) XI^#CSTA^ 
format (2F10.A) 

NIN . ABS(XIN) 

CALL GOPLCT(NRN) 

IF I MN.GT.NINMaX ) GU TO 90 
READ IN 0(S) 

OMAX » 0. 

DO 20 J « 1/NlN 
READ(N6»510) SF(J)#0I(J) 

F0RMAr(2E20.10) 

QMAX a AMAXl(QMAX»ABS<Ol(J)n 
CONTINUE 

XFAC » 2./(SF(NlN»-SF(l)l 
CONST a -l.-XFAC*SF(l) 

FAC a 1, 

I F ( QM AX • GT • 1 • 6 ) FAC * l«o/QMAX 
WFIFE (N<t,130) 

DO 3C J a 1/NlN 

SX(J) a XFAC*SF(J)+CCNbT 

QX(J) a FAC*01(J) 

= XMACH(QX(J)*uX(J) ) 

Cl( J^SF(J),0I( J)^SX( J),QX(JJ,0{J) 

CONTINUE 

WRITE (N'f»oOC) CSTAR 

SIZE a .lA 

sex a 5. 

SCY a 2.5 
XGR a 6,0 

IF (XIN.LT.O.) XOR a p.7p 
CALL PLUT( ADR, 5, ) 

DO 60 Lal#NIN 

PL0T(-5.0»6.0»2) ’ 

SYMBOL(-<».5#3.5»SIZE/lHO,0.il» 

SYMBOL (-5.0»SCY*CSrAR,2.*SI2E> 15#0.,-1) 

S YM30L ( -5 .0»-SCY*CSTaK» 2. *SIZE» 15»0.,-I ) 

L = l>9 


CALL 
'■ALL 
CALL 
CALL 
Call 
DO 7C 


70 


Yh a FL0AT(L-1)-A.0 

S a ,4*FLUAT(L-1)-1.6 

CALL SYMBUL(-5.C>YH, SIZE>15,0..-1) 

ENC0u6(10,100/A) S 

CALL SYMbOL(-5.7* Yri,SIZF,A#0.#<t) 


.)*QS)) 
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eo 

100 


‘tO 


95 


CALL PLOT( “5 • 0> 0»> 3 ) 

CALL PL0T(5.0»0.»2) 

XH » FL0AT(L-1)-5.0 
S « .2<-FLCAT( L-D-l.O 

^'’”0 uisiribution,o.,2a) 

DS » ( SX(NIN) -SX(1 ) )/FLOAT(NCPT-l) 

FIND Q(S) AT EVENLY SPACED POINTS 
ES(1) = SX(1) 

DO AC J = 2fNQPT 
ES(J) « ES(J-1)*0S 
ES(NCPT) « SX(NIN) 

CALL SPLIFCNIN # SX, OX, GP, GPP, GP PP, 3» C. , 3, C . ) 

CALL INTPL(N4pr,ES,0,SX,QX>GP,GPP»6PPP) 

CALL PL3T(SCx*ES(1),SCY*Q(1),3) 

DO 95 L=2,N0PT 

CALL PL0T(SCX*bS(LJ,SCY*Q{L),2) 

CALL PLOT (-X0R>-5. 5,~3) 

CALL FRAME 

IF (CSTAR.LT.O. ) 60 TO 210 

CALL SPLIFCNCPf, es,0,6P,GPP,6PPP,3>0.,3,0.) 

CALL INTPLI(1,SC, 0«,NOPT>£S,0,GP,GPP/6FPP) 

INTEGRATE Q(S) TO GET PHI(S) 

CALL SPLIE(NQPT,ES,Q,O00S,0PP,PHT,-3,0.,3,0.) 

Call intpl(i>so, phhn,es,pht,q,dods,opp) 

GAM « PHT(NOPT)-PHTd) 

SCALOI » PHT{ NQFT )-PHMN 
DO 50 I » 1,N0PT 
OKI) = PHT(I)-PHMN 
VAL = AMAXI (0 Cl ( 1 ) ) /SC ALOI 
PHT(I) . SI6N(S0PT(VAU,ES(I)-S0) 

REWIND N6 

RETURN 

WRITE {NA,110) NINHAX 
CA<,L PL0r(C.,0.,999) 

CALL EXIT 
RETURN 

*^^'^‘^'2F15.6,3X,2F15.6,3X»F15.6) 

END 


50 


90 

210 


110 
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SUBROUTINE INCOHP 

« /RI>M31),Sfi01, ““’"‘'^“'*“'"'<C0l“n62),0tl.01.0(lolr 


10 


2 't'"C5.C6,C7,B£;;;[lA;fstj;j;;i'f‘;J5“;J«;A2SU,0CRlI,cl.CA 

3 > JK» I2» ITYP/MODE# IS»NFC/NCY N«n ' ^ ^ ^ NiP 

^ *'^PIS,LLW,LSCP,h;;Uw;£PS 1 

i ^SCALOD,hc,GA«HA,NOFr.CST*i 

'j;?C,QFAC,GAn,KO£S,kTsLo%:'o?i 

OFAC . SCAL0I-.5*GAM 
DO 1C 1=1,10 

TaU = ASlN(-GAM/{Pl*Qf AC) ) 

TAU*^. ASINfiGAM/fpiIoFAClT'^^^^^^^'^^^^^*^ 

BBd) » -TAU-ALP 
DPHI s A.*GAN/0FAC 
CALL COSI 
ANG = 0. 

DO 20 I-1,M 

PHI(KH,1) . PHI(1,1)+^,aM 

PHI(PM+1,1) « PHI(2,1)+GAM 

INC = 1 

CALL CYCLE 

PETUFN 

END 


bUBPOUTINE CYCLE 

rf Jm? y » i'^DC 0 

2 >SI( 162 ),PHIP( 162 ),xc l 62 LYrMljy:y”y’'®®‘^«> 2 ),CC{l 62 ) 

2 'C 4 ,C 5 ,C 6 ,C 7 ,bET,bElA,FSYM; 5 sE? y?M 

3 > IK, JK, d, ITYP,MCDI, IS,NFC»tNCY NRM I 1 ^ ^ ^ ^ NSP 

^ ^ N 6 # o A h^*A > NOP C S T Aft • P E- M nco ^ t i* ^ 

6 ,INC, 0 FAC, 0 An,KDE 3 ,PUsL"pu^?;‘’"’“‘''"’'"'"'’‘‘'''' 

^ '“P'i'>Am),A 2 (ii,A 3 n;,«}u"j;i;;’’“'*’'““''SA.i).Bsm 
1 
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? ‘P'*PP(1>>PP(1»15))/(FPPPP(1),fP(1, 17)) 

3 'J®^^^*'PPtl>19))»(PHIv(l)>FP(l/21))i<PHlS(l)»FPn.?'-.»i 

A *1 *^^^^^*^**' *Am)»f^P(l) )> (A2(1)»RP(7) ) 

data KO/0/ 

LC ■ NFC 
NMP « 2*tc 
HC « NMP ♦! 

PILC « PI/FLaAT(LC) 

NCUm » NHP/M 
NDUM2 • NLUHl-1 
IF (INC.EO.l) GO ;o 60d0 
INOCD « 1 

CALL GTURe( 0 .» 0 .» *4»CDx»0.».125>,l) 

INOCO s 0 
6080 DO 10 I»1>MH 

CIRCd) a FLOAT(I-l)*tJT 
10 PHIS(I) * PHI(I#1)+C0(I) 

Bbl « bS(l) 

00 6 I»1,MM 
0(1) a PHIS(I) 

6 SS(I) a CIRCd) 

DO 7 I«1,HC 

7 CIRC(I) a FLOAT( I-1)*P1(.C 

CALL INTPL(MC»C1RC#PHIS>SS>U/CPHDW» FPPP»FPPPP) 

5 00^9 FPPP, FP»PP,1/0.,1,0.) 

SS(I-tO) a CIRC(I) 

9 Q(I-40) a OPHOW(I) 

CALL SPLIF ‘80»SS»Q»FPP>FPPP,FPPPP>3»0.,3»0. ) 

CALL lNTPLI{l#wNP^O#/dOfiS>0/FPP^FPPP^FPPPp) 

if‘^*C^PC»PHIS#DPHDWAC2PHDW>FPPP,l» 0.>l,0.) 
cr ‘*^'*P» PMP*N#CI(<C>PHIS> OPH1)W>02PHOW> FPPP) 

SCALCO a PH1S(HC)-PHHN 
REWIND No 

READ (N6) (PHlTd),SS(n,PHlV(l),0(I)>I-l,NUPT) 

(jALL SPLIMNQPT>PHIT>Oj>FPP>FPPP>rPPPP>3*0.>3#0.) 

VAL a SORT((PH1S(1)-PHMN)/SCALQO) 

DO 20 I«1,MC 

VAL a AMAX1(0.,PHIS(1)-PH.1N)/SCALQ0 
FAC a 1, 

IF (CIRC(I).LE.WNP) FAC « -1. 

20 PHIS(I) a FAC*SQRT(VAL) 

CALL INTPL(HC»PhIS,Qx>PHIT>C» FPP» FPPP.FPPPP) 

IF (INC.EO.l) GO TO 6837 

determine Em 

XNUM a 0. 

DEN a 0« 

DO ^A4A J a 2/M 


)♦()$)) 


160 







444b 


K « J 

IF (MC.NE.MM) KaNDUMl*J-NDUM2 

FPPP(J) . VAL 

> CO TO 4^4^ 

XNUM = XNOM + QX(K)*CX(K)*VAL 
DEM « OEM + VAL»VAL 
CONTINUE 

OINF » SuRT(XNUM/DEN) 

DQMAX » 0, 

DOAVt a 0. 

DO 444b J a 2#M 
K a J 

IF (HC.NE.i'IM) K = NDUN1* J-N0UM2 

C7 ■ GAMMA/(6AMMA-1 . ) 

EMX a EM 


EM 

LM 

EM 

EM 

Cl 

C6 

C4 

Cb 


(GAMMAn.)*CSTARtCSTAR/(OINF*OINF)-GAMMAn. 

SORKEM) 

{ 1 •”PEM) *EM+REM*tMX 
CZ+l./IEM+EM) 

Ca^EM+EM 
1 • ♦Cb 

l./(C6*C7) 

QCRIT = (Cl+Cl)/(GArtMA+l . ) 

determine scaling for phi by least SOUARFC pit 

XNUM = 0. 

DEN . 0. 

00 ^7 
K a J 

IF Cf^C*Nt.MM) Ks NDUMl ♦ J-,^0UM2 
XNUM a ANUM+PHI(J,l)+Cd(j)-PhMN 

DEN a DEN+ Q(K) 

FAC a XNUM/OEN 
OPHI a FAC»GAM 
8887 DO 0d86 J*1»HM 

K a J 

flftflA KaNDUMl#J-NUL'M2 

8886 CCP(J) = OX(K) 

CALL SPLIF(MOPT»PHlT»SS»FPP»rPPPi FPPPP » n o « . 


47 
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^0 CONTINUE 

VAL « A3SIPPPP(1)) 

OSilJ ■ ALOGIVALI 
US(MC) » 05(1} 

A3 CONTINUE 

DD 2A0 1 -1,HC 

ANGL ■ FLOAT! I-l )*pxLC 

CXdl a COS(ANGL) 

2A0 SXd) a SIN(ANGL) 

CD XC40 Ial,LC 
FPPP(I» a AA(I) 

FPPPP(l) a BB(I) 

AA(I) = CX(2»I-1) 

1040 fcB(I) a -iX(2*I-l) 

CALL F0UCF(NHP#LS>FPP>AA>8B) 
DFAVE a 0. 




= 1 >LU 


AA( I ) = -aA( I) 

IF (FPPP(l).EO. 9999 V.) GO TO lu 50 

Sb!!,’ : 

DFAVfc a SQRT(0FAV£/160. ) 

BBd) a Bbl 

IF (INC. to. 0 ) GG TO 45 
QINF a OFAC/(^.*tXP(AAd)) ) 

.5*(GAK,'*A-1 . ) 

gamma/ CGaNNA-1.) 

Cl a C 2 + 1 ,/(FM*EM) 

C 6 a C 2 *EM*EM 
C4 a 1 . tC6 
C 5 a 1 ,/{C 6 *C 7 ) 

WCRIT a (C1+C1)/(GAHMA+1.) 

A5 CONTINUE 

IF (INC.tC.l) GO TO 6030 
IF (KDES.EO.l) g:) TO o040 

tOPO ”lt£"?SL6S:S?“ ’•h£.JI.ANC.(KO.Ut. 2.) (.0 TQ tcCL 

toco ANGO a -PAO*Bdd) 

AL°1 a ALP+PAO 


C2 

C7 


1 rC * *^C» 0 UAOt,D 3 HAXf LFA'/ 

f 020 PORM/T ( IX, 13,5(2x,fQ.j),2;;,I3,^X, 

1 » 2X, F3 . 5, 2X, F5. 3 ) 

6030 INC a 0 


Vk, r A, N5P, tM, CL, ,.LP1, A NOC 
• 4,3 {aiX,F7.4) 


ir- 
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KD « KO+1 

IF (KL.viT.NOfcS) KC»1 
CALL MAP 


fOG 


805 


; 5 t >5 


CL = 2.*uPHI*CH0 
KtWIKO ilA 

I fb ( M<i# 0 00 ) MM 
FCRMAT(10X,I3) 

SNX = 1, 

*<RITb TC»uJAVb 

Format (F 7 . 3 , 2 tio. 2 ,F.,n 

JG 5555 I ■ 1,MH 

VAL » XM4CH(CCP(I)«CCP(I)J 

CCP(l) » CP(VAL) 

wRITb(M<»,t)10) XC(I),CCP{I) 

CuNTINUE 


blC FCR1AT(2F10.A) 
CALL COSI 
tPSl » EP51-0EP 
RETUFN 
END 


# YR,SNX 


2 

3 

A 

5 

6 


100 


120 


iuaKu'UTlNE OOTPT 

ii:?;5SSSs:sS[;£;^^ 

PE WIND N3 
Xf: a M*< 

IF ( X ju r . (j r, 0 . ) gl to ijo 
00= SORKCCPIT) 

WP I TE ( N3 » 1 20 ) XM»OQ 
F oRMAT ( 2P1v)*A ) 

DO lAO L= 2 ,M 

0 = (PHnL*l,U-PHl(L-i,l) )ntL iM-bKU 

OS = (0*0) /F->( L, 1) 

C'(L> = SOPKJSl 
?HIS{L- 1 ) = Phl(L,n*CU(L) 

ClRC(L-l) = FLuAT (L- 1 J»Or 


lAC CuNTlNUE 
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'iiSKar 


, Tin: 

IS POilit 



0(1) » .5*{Q(2)+0(n)) 

UtHM) * 0(1) 

HZ « H-1 

CALL SPLIMMZ»CUC»PHIi,FPP»FPPP, FPPFP, 3»0 . > 3/ 0. ) 
CALL S«»LIF(MZ#ClRC*hPP»FPPP/FPPPP,FPPPP?»3#0.#3»0.) 
CALL INTPLI(l»WNP>0.>f1Z»CIRC»FPP>FPPP#fPPFP,FPPPPP) 
0 ( 1 ) » - 0 ( 1 ) 

DO 200 1=2, M 

0(1) = SIGN(0(1),CIRC(I-1)-WNP) 

200 CONTINUE 
C Ci(MH) « 0(MH) 

ARCl » 0. 

WRITE(N3,160) AkCl,0(l) 

PRINT 160, ARCl, 0(1) 

00 160 L=2,HH 
DX « XC(L)-XC (L-1) 

CY » YC(L)-YC(L-1) 

APCl = Af<Cl ♦ SORT(DX*OX + OY*CY) 

WPITt(N3,lbO) ARC1,Q(L) 

PRINT 160, ARC1,0(L) 

160 F0RHAT(2E20.10) 

150 CONTINUE 
XDUT « 0. 

RETURN 

END 


0 


BLOCK DATA 

COMMON PHI(162,31),h?(162,31),A(31),B(31),C(31),D(31),E(il) 

1 »PP(31),PPP(31),R(31)»RS(31),l<I(3i),AA(lo2),u3(lc2),CL(io2) 

2 »SI(162),PHIR(162);XC(162),YC(lb2),FM(ltj2),AkCL(162),L'SuM{lL2) 

3 #ANG0L0(162) ,XUL0(162), Y0LD{162), ARC0LJ(lt.2),0EL0L0(lt2) 

,RP4(3i),RP5(31) 

COMMON /A/ PI»TP,RAO,tM,ALP,RN,PCri,XP,TC,CHL»,uPHi,CL,hCL,YR 

1 »XA,YA»TE,OT,OR,CELTri, CELR jFA, LCN,OSN,RAA,tPSIL,OCRlT,tl»CC 

2 /CA,C5,C6,C7,3ET,bETA,FS YM,XSLP, SbPM, TT lE (A ), M,N,MM, NN,N6P 

3 , IK, JK, IZ, ITYP,MOl)E, I NFC, NC Y, NRN, NG, I J I M, N2, N3, NA, NT, 1 X X 
A ,NPTS,LL,I,LSEP,MA,NE/i(, bPSl,NL>ES,XLtN,iCALQl 

5 ,SCALQ0,N6,GAHMA,N0PT,CSTAR, PLM, OLP, wINF, T iTtP, XOCT 

6 ,INC,OFAC,GAM,KL£S,PLTS2,OPL,OPO 

•♦♦♦IDIM MUST BE SLT TJ THE FIRST OIELNSILK uF Phi**** 

DATA PI/3.iA169265i5c97-)/ , tM/.7b/ , ALP/0./ , CL/lCO./ , 

1 PCH/,07/ , FSYM/1.0/ , RCL/1.0/ , otTA/0.0/, RN/2C.cc/ » 

2 S6PM/.00A/ , XSI.P/.9J/ , <P/O.C/ , M/lbO/ , N/30/ , ^^N/l/ , 

3 NFC/80/ , NPTS/81/ , LL/0/ , NG/1/ , li/2/ , I0ih/lt2/ , MODE /!/ 
A , JK/0/ , N2/2/ , N3/3/ , NA/A/ , LSEP/lbi/ , 1//126/ , IlYF/1/ 

5 ,N0PT/321/,REH/0./»kPil/0./,CS1AP/lC0./ 
t> ,DbP/0./,N0ES/-l/, TsrtP/. 2 /»KDLS/x 0 /,FLTSt/» 0 ./, 0 PL/.t> 5 /, 

7 CPU/. 96/ 

END 
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